% GCI.m - example how to calculate the Grid Convergence Index (GCI) [1]
clear, clc;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% INPUT SECTION
% Number of cells for different grids, chosen such that (r21 AND r32) > 1.3
N1 = 500000;
N2 = 150000;
N3 = 50000;

% Values of property for different grids, phi
phi1 = 6.063;
phi2 = 5.972;
phi3 = 5.863;

% Total volume of the computation domain for different grids, V [m^3]
V1 = 1E-2;
V2 = 1E-2;
V3 = 1E-2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Representative grid size for different grids, h
h1 = (1/N1*V1)^(1/3);
h2 = (1/N2*V2)^(1/3);
h3 = (1/N3*V3)^(1/3);

% Refinement factor r
r21 = h2/h1;
r32 = h3/h2;

eps32 = phi3 - phi2;
eps21 = phi2 - phi1;

% Apparent order, p 
% p = 1/log(r21)*abs(log(abs(eps32/eps21)) + ...
%     log( (r21^p-sign(eps32/eps21)) / (r32^p-sign(eps32/eps21))));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Fixed point iteration method
g=@(x) 1/log(r21)*abs(log(abs(eps32/eps21)) + ...
    log( (r21^x-sign(eps32/eps21)) / (r32^x-sign(eps32/eps21))));

% Start out iteration loop
x1 = 1/log(r21)*abs(log(abs(eps32/eps21)));
x2 = g(x1);
it = 0;     % iteration counter
while (abs(x2-x1) > 1e-5 && it<100)
    it = it + 1;
    x1 = x2;
    x2 = g(x1);
end
% fprintf('#iterations: %d\n', it);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
p = x2;

phi21_ext = (r21^p * phi1 - phi2)/(r21^p - 1);
phi32_ext = (r32^p * phi2 - phi3)/(r32^p - 1);

% Approximate relative error, e21_a
e21_a = abs((phi1 - phi2) / phi1);

% Extrapolated relative error, e21_ext
e21_ext = abs((phi21_ext - phi1) / phi21_ext);

% Fine grid convergence index, GCI21_fine
GCI21_fine = 1.25 * e21_a /(r21^p - 1);

pad = '--------------------';
fprintf('%s\n',pad);
fprintf('Input in SI-units: \n');
disp(table(N1, N2, N3, V1, V2, V3, phi1, phi2, phi3));
% disp(table(phi1, phi2, phi3));
fprintf('%s\n',pad);
fprintf('Output to check: \n');
disp(table(h1, h2, h3, r21, r32));
fprintf('%s\n',pad);
fprintf('Report the following: \n');
disp(table(p, e21_a, e21_ext, GCI21_fine));
fprintf('%s\n',pad);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% References
% [1] I. B. Celik et al., J. Fluids Eng. 130, no. 7 (2008)