[esha]

Hi, I want to calculate the GCI. I have three different mesh sizes (M1, M2 and M3 with 7227990, 9133110 and 9761535) elements.
Mesh Cd Cl
722799 0.1651 0.4184
913311 0.1662 0.4180
976153 0.1661 0.4180
y+ =0.02 and spacing of first cell is 0.0003.
I can calculate the error by using formula error =(f2-f1)/f1. But unable to figure out how to calculate GCI. Can anyone please explain and help in this regrad? Thank you!

[LuckyTran]

The NASA page has an excellent example, although not a thorough as Roache's original text.


You need the grid refinement ratio r and order of convergence p and then

[Ford Prefect]

I have always found it a bit difficult to interpret which type of estimation is used when GCI is reported in various texts.

GCI_23 is sometimes also reported as GCI_32, while still estimating on the “fine” grid (i.e. grid level 2). Similarly with GCI_12 and GCI_21. The second point being that having three (or more) grid levels and reporting a GCI_23 where level 2 is the “fine” level is a bit confusing. Perhaps the terminology of “fine” and “coarse” should be dropped entirely as that was more relevant in the original text where only two grid levels were used?

[LuckyTran]

The original text by Roache uses three grids...

GCI 21 or 12 as well as 32 and 23 refers to whether you are estimating the error by extrapolating from the coarser grid or from a finer grid. There is no rule that mesh sensitivity studies must always be performed by using a finer and even finer grid. One can do a coarser and even coarse grid or one coarse and one finer. All 6 GCI's are therefore valid.

Furthermore, there is no rule that there can only be 3 grids. It is even better to use 4, 5, 6, N grids.


1+2=3 is confusing for some people that can't differentiate between a + and ÷ symbol. Or they can just learn the difference.

[Ford Prefect]

Originally Posted by LuckyTran

GCI 21 or 12 as well as 32 and 23 refers to whether you are estimating the error by extrapolating from the coarser grid or from a finer grid.

Yes, and some report GCI_21 and use the “fine” estimation. For instance, in the following reasonably recent text by Roache and others [DOI:10.1115/1.2960953]

Here they report GCI_fine^21, would you say that this is consistent with the NASA text or your explanation above?

Originally Posted by LuckyTran

There is no rule that mesh sensitivity studies must always be performed by using a finer and even finer grid. One can do a coarser and even coarse grid or one coarse and one finer. All 6 GCI's are therefore valid.

Furthermore, there is no rule that there can only be 3 grids. It is even better to use 4, 5, 6, N grids.

1+2=3 is confusing for some people that can't differentiate between a + and ÷ symbol. Or they can just learn the difference.

Thank you for these insights. I agree.

[LuckyTran]

1. If you are doing a one-off simulation, then you often refine the grid and present the results from your finest mesh, because that is the one that is most likely to be most accurate. "omg look at this tiny wall y+ value that I get!"
2. If you do a DOE study, you estimate your error once and then run a thousand cases using a coarse grid because you can't afford to a super fine grid every time. Or maybe you can, that would be exceptional.

Both scenarios are explicitly mentioned in Roache's text and the NASA blog. Now, the vast majority of CFD publications are from the 1st category and so there is a general preference for reporting the result one way (if at all).

If I give you two consecutive terms in a sequence, say: 1 and 4. I can ask you, what comes after 4? And what is the error in your estimate of the term after 4. Or I can ask you what came before 1? And what is the error in the estimate of that term? That is GCIfine and GCIcoarse

[Ford Prefect]

Originally Posted by LuckyTran

1. If you are doing a one-off simulation, then you often refine the grid and present the results from your finest mesh, because that is the one that is most likely to be most accurate. "omg look at this tiny wall y+ value that I get!"
2. If you do a DOE study, you estimate your error once and then run a thousand cases using a coarse grid because you can't afford to a super fine grid every time. Or maybe you can, that would be exceptional.

Both scenarios are explicitly mentioned in Roache's text and the NASA blog. Now, the vast majority of CFD publications are from the 1st category and so there is a general preference for reporting the result one way (if at all).

If I give you two consecutive terms in a sequence, say: 1 and 4. I can ask you, what comes after 4? And what is the error in your estimate of the term after 4. Or I can ask you what came before 1? And what is the error in the estimate of that term? That is GCIfine and GCIcoarse

Ok, so how would you report the GCI on the medium grid in a sequence of 3 (the h2 grid, with h1 being the finest grid and h3 being the coarsest)? Can you think of two ways of doing this? If yes; is it good if such reports are consistent so that the reader understands which type of estimate was actually used without having to redo the calculations? I am not asking for help how to interpret the method. I just stated that there is a confusion in the literature on how to report values, with regards to the fine and coarse estimates.

We have, in different texts

GCI_12
GCI_21
GCI_fine^21

All with the same meaning, using the fine estimate. The GCI_fine^21 is perhaps the best one since it explicitly states that the fine estimate is used. The order 12 or 21 is irrelevant if “fine” is also added as subscript.

Any thoughts on this?

[LuckyTran]

You are arguing semantics on whether a medium-well cooked steak is medium-well or well-medium, more-than-medium, or less-than-well. There is no need to be dogmatic. Not everyone uses the same English, and (statistically) most people don't even English. There is a reason a nomenclature section is required.

The real question is... when does it even matter? 99% of CFD is presented without error estimates in the first place.

GCI is rooted in statistics. Even without employing Richardson extrapolation, there are at least 3! factorial ways of estimating the error. Be glad that Roache has given a prescription to limit this to only two...

Even mathematicians can't agree on how to define concavity...

[Ford Prefect]

Originally Posted by LuckyTran

You are arguing semantics on whether a medium-well cooked steak is medium-well or well-medium, more-than-medium, or less-than-well. There is no need to be dogmatic. Not everyone uses the same English, and (statistically) most people don't even English. There is a reason a nomenclature section is required.

The real question is... when does it even matter? 99% of CFD is presented without error estimates in the first place.

GCI is rooted in statistics. Even without employing Richardson extrapolation, there are at least 3! factorial ways of estimating the error. Be glad that Roache has given a prescription to limit this to only two...

Even mathematicians can't agree on how to define concavity...

I agree that a nomenclature section is very useful. If your field of CFD ignores error estimates in 99% of the cases then I would argue that there is room for improvement. I would never accept any CFD paper without any form of grid sensitivity analysis. I would never accept any experimental paper, that is based on a one-time experiment either.

I do like the effort by Roache. I just pointed out that there is a problem with consistency for the (1%?) of papers that do an analysis and I wished to have a simple discussion about it. I get the feeling that your passive aggressive response indicates you do not wish to discuss this matter

[LuckyTran]

And you are the editor of which journal? How many papers have you rejected?

It is a simple fact that most papers don't present error analysis. And even when they do, it is never clear whether they use a inverse z score or an inverse t score or whether they correctly use the sample statistic versus the population statistic. You are hijacking someone else's post about how to calculate a GCI because you have your own agenda to push because you get too easily confused by two GCI's where there are an infinitude of other potential definitions for grid error estimates. If you have a problem with the way people present error analysis, go tell that to the journal editors! Anyone can nitpick behind a keyboard

[Ford Prefect]

No need to turn this thread into a common toxic Reddit thread. If you do not wish to discuss this topic, then nothing is forcing you. Just move along.

I wrote up a small code in Julia for anyone that has not been deterred by the level of the thread so far.

# Function to calculate GCI
# References:
# DOI: 10.1115/1.2960953
# https://www.grc.nasa.gov/www/wind/valid/tutorial/spatconv.html

function calculate_GCI(N, ϕ, ΔV, dimensions)
    # N[]: number of cells at grid level i, where N[1] is the finest and N[3] is the coarsest, select exactly three grid levels
    # ϕ[]: Some value of interest at grid level i
    # ΔV: Domain volume (or area for 2D)

    # Array consistency check
    if (length(N) != 3) || (length(ϕ) != 3)
        println("Wrong array length.")
        return
    end

    # h is based on the simple average cell volume (or area for 2D)
    # Step 1. Define a representative cell, mesh, or grid size h. For example, for three-dimensional calculations,
    if dimensions == 3
        h = (ΔV ./ N).^(1/3)
    elseif dimensions == 2
        h = (ΔV ./ N).^(1/2)
    else
        println("Wrong dimensional input")
        return
    end

    # Select three significantly different sets of grids and run simulations to determine the values
    # of key variables important to the objective of the simulation study, supplied here as the array ϕ

    # Calculate grid refinement factor, it is desirable that it is higher than 1.3 (from first reference)
    r_21 = h[2] / h[1]
    r_32 = h[3] / h[2]

    # -------------------------------
    # Calculate the apparent order, p

    ϵ_21 = ϕ[2] - ϕ[1]
    ϵ_32 = ϕ[3] - ϕ[2]

    p = 1 / log(r_21) * abs(log(abs(ϵ_32 / ϵ_21)))

    residual = 1
    while residual > 1e-6
        s = sign(ϵ_32 / ϵ_21)
        q = log((r_21^p - s) / (r_32^p - s))
        p_new = 1 / log(r_21) * (abs(log(abs(ϵ_32 / ϵ_21))) + q)
        residual = abs(p_new - p)
        p = p_new
    end

    # -------------------------------


    # Extrapolated values
    ϕ_ext21 = (r_21^p * ϕ[1] - ϕ[2]) / (r_21^p - 1)

    # Error estimates
    ϵ_a21 = abs((ϕ[1] - ϕ[2]) / ϕ[1]) # approximate relative error
    ϵ_a32 = abs((ϕ[2] - ϕ[3]) / ϕ[2]) # approximate relative error
    ϵ_ext21 = abs((ϕ_ext21 - ϕ[1]) / ϕ_ext21) # extrapolated relative error
    
    GCIfine_21 = 1.25 * ϵ_a21 / (r_21^p - 1)# how far grid 1 is from the asymptotic value.
    GCIcoarse_21 = 1.25 * ϵ_a21 * r_21^p / (r_21^p - 1) # how far grid 2 is from the asymptotic value.

    GCIfine_32 = 1.25 * ϵ_a32 / (r_32^p - 1) # how far grid 2 is from the asymptotic value.
    GCIcoarse_32 = 1.25 * ϵ_a32 * r_32^p / (r_32^p - 1) # how far grid 3 is from the asymptotic value.

    return (r_21, r_32, h, p, ϕ_ext21, ϵ_a21, ϵ_ext21, GCIfine_21, GCIfine_32, GCIcoarse_21, GCIcoarse_32)
end

# Example input arrays from reference DOI: 10.1115/1.2960953
N = [18000, 8000, 4500]
ϕ = [6.063, 5.972, 5.863]
ΔV = 1 # does only matter if the domain volume actually changes with changing grid size

# Calculate GCI
(r_21, r_32, h, p, ϕ_ext21, ϵ_a21, ϵ_ext21, GCIfine_21, GCIfine_32, GCIcoarse_21, GCIcoarse_32) = calculate_GCI(N, ϕ, ΔV, 2) # note that reference DOI: 10.1115/1.2960953, is a 2D case

# Simple output
println("r_21: $r_21")
println("r_32: $r_32")
println("Representative cell sizes (h): $(round.(h, sigdigits=3))")
println("Apparent order (p): $(round(p, sigdigits=3))")
println("Extrapolated value (ϕ_ext21): $(round(ϕ_ext21, sigdigits=3))")
println("Approximate relative error (ϵ_a21): $(round(ϵ_a21*100, sigdigits=3))%")
println("Extrapolated relative error (ϵ_ext21): $(round(ϵ_ext21*100, sigdigits=3))%")
println("Grid Convergence Index (GCIfine_21): $(round(GCIfine_21*100, sigdigits=3))%")
println("Grid Convergence Index (GCIcoarse_21): $(round(GCIcoarse_21*100, sigdigits=3))%")
println("Grid Convergence Index (GCIfine_32): $(round(GCIfine_32*100, sigdigits=3))%")
println("Grid Convergence Index (GCIcoarse_32): $(round(GCIcoarse_32*100, sigdigits=3))%")

[LuckyTran]

There is already a fortran code provided on the decades old NASA thread

You go beyond the level of micromanagement and write code to force people to do GCI your way like a true tyrant and you call me the toxic one!? Oh please. Just because you couldn't convince some redditors to use your nomenclature... I offer ppl freedom to choose the GCI as they like

[esha]

Thank you so much guys for your help.

[Ford Prefect]

There is nothing wrong using the Fortran code from NASA. It is, however, assuming constant grid refinement ratio. In case you are using unstructured and somewhat automated meshing where it is impractical to enforce constant grid refinement ratios then it is better to use the method suggested by DOI: 10.1115/1.2960953, this reference is written by Celik, Ghia, Roache, Freitas, Coleman and Raad. I believe this is a good reference.

The Julia code I provided is just based on that reference. I have not added any "personal nomenclature", as I do not have any such, which I thought was obvious from my inputs earlier. All nomenclature is from the reference. The only thing I have added is the "coarse" values, as calculated by the NASA (Roache) reference.

[esha]

Originally Posted by Ford Prefect

There is nothing wrong using the Fortran code from NASA. It is, however, assuming constant grid refinement ratio. In case you are using unstructured and somewhat automated meshing where it is impractical to enforce constant grid refinement ratios then it is better to use the method suggested by DOI: 10.1115/1.2960953, this reference is written by Celik, Ghia, Roache, Freitas, Coleman and Raad. I believe this is a good reference.

The Julia code I provided is just based on that reference. I have not added any "personal nomenclature", as I do not have any such, which I thought was obvious from my inputs earlier. All nomenclature is from the reference. The only thing I have added is the "coarse" values, as calculated by the NASA (Roache) reference.

Yes, I have gone through the page of NASA, and they have a constant spacing ratio 2. Thank you for suggestion. I will read the references( you mentioned), and will get back soon. Thank you once again. If will have some issue about GCI, I will ask in this thread again.