% Lid driven flow convergence test
clear all, clc, clf, close all 
%Set number of timesteps
app.dt = 1/400;
%Set plot intervals
app.PlotIntrvl = 4000;
%Set viscosity 
app.nu = 1/1000; 0.005; %(for RE=200)
U=cell(4,2);
for i = 1:length(U)
    N = 4*2^i; 
    app.Nx = N;
    app.Ny = N;
    [P, u, v, rho] = TestCase(10,3,app);
    %extract center u and v velocities
    if ~mod(N,2)
        U{i,1} = u(5+N*(N/2-1):4+N*(N/2));
    else
        U{i,1} = 0.5*(u(5+N*(N-1)/2:4+N*(N-1)/2+N)+...
            u(5+N*(N+1)/2:4+N*(N+1)/2+N));
    end
    dy = 1/N;
    U{i,2} = fliplr(linspace(dy/2,1-dy/2,N));
end
%% Plot figure
figure
str={'--','xg','or','<m','*c','k'};
hold on
for i = 1:length(U)
    hold on
    plot(U{i,1},U{i,2},str{i})
    ylabel('y')
    xlabel('U - velocity')
    title('Vertical center line U-velocity')
end
legend('8x8','16x16','32x32','64x64','128x128','256x256')
%% Convergence
%Compare everything on fine grid
x = U{end,2}';
for i = 1:length(U)-1;
   err(i) = norm((U{end,1}'-interp1(U{i,2}',U{i,1},x','spline','extrap')))./sqrt(length(x)); 
end
log(err(1:end-1)./err(2:end))/log(2)
figure
semilogy(err)
hold on
semilogy([err(1), err(1)/4, err(1)/16], 'g')

