% clear all, clc, clf, close all
function [] = LidDrivenCavityConvergeScript_num(Re,runref)
if nargin < 1 || isempty(Re)
    Re = 500;
end
if nargin < 2 || isempty(runref)
    runref = true;
end

%Setup generics
d0=fileparts([pwd,filesep]);
d0=d0(1:strfind(d0,'trunk')+4);
sscript = 'LidDrivenCavitySetup_num';
pscript = 'LidDrivenCavity_num_Plots';

%Create the app for the Reference solution
pscriptref = 'LidDrivenCavity_num_refPlots';
app.DO_cor_iter = 0;
app.Nx = 512;
app.Ny = app.Nx;
NNXX = app.Nx;
NNYY = app.Ny;
app.dt = 1 / app.Nx / 8;
app.T = 5;
app.PlotIntrvl = round(32 * app.Nx / 128);
app.nu = 1/Re;
app.S = 3;
app.savemat = ['DOnumPap/LidRefHR',num2str(Re),'_',num2str(1/app.dt)];
%app.saveimg = app.savemat;
app.MC = 10000;
app.var = 1;

%Run the reference solution
if runref
    NavierStokesInc_DO_mx(app,sscript,pscriptref);
end
beep
%Run the lower-res cases
if isfield(app,'saveimg')
    app = rmfield(app,'saveimg');
end
app.ref = app.savemat;
app.dt = app.dt / 2;
app.PlotIntrvl = app.PlotIntrvl * 2;
for j = 1:4
    app.dt = app.dt * 2;
    app.PlotIntrvl = app.PlotIntrvl / 2;
    app.savemat = ['DOnumPap/Lid',num2str(Re),'_',num2str(1/app.dt)];
    for i = 1:5
        app.Nx = NNXX / 2^(i - 1);
        app.Ny = NNYY / 2^(i - 1);
        %%
%         app
%         NavierStokesInc_DO_mx(app,sscript,pscript);
        beep
        if exist(sprintf('%s/Save/%s/RE%d_NX%03d_dt%03d.mat',d0,app.savemat,1/app.nu,app.Nx,1/app.dt), 'file')
            load(sprintf('%s/Save/%s/RE%d_NX%03d_dt%03d.mat',d0,app.savemat,1/app.nu,app.Nx,1/app.dt));
            fprintf('%s/Save/%s/RE%d_NX%03d_dt%03d.mat\n',d0,app.savemat,1/app.nu,app.Nx,1/app.dt)
            %Calculate final norms
            PL2f(i,j) = sqrt(sum(PL2) * app.dt * app.PlotIntrvl); %#ok<SAGROW>
            PiL2f(:,i,j) = sqrt(sum(PiL2,2) * app.dt * app.PlotIntrvl); %#ok<SAGROW>
            uL2f(i,j) = sqrt(sum(uL2) * app.dt * app.PlotIntrvl); %#ok<SAGROW>
            uiL2f(:,i,j) = sqrt(sum(uiL2,2) * app.dt * app.PlotIntrvl); %#ok<SAGROW>
            vL2f(i,j) = sqrt(sum(vL2) * app.dt * app.PlotIntrvl); %#ok<SAGROW>
            viL2f(:,i,j) = sqrt(sum(viL2,2) * app.dt * app.PlotIntrvl); %#ok<SAGROW>
            rhoL2f(i,j) = sqrt(sum(rhoL2) * app.dt * app.PlotIntrvl); %#ok<SAGROW>
            MEAN_ENERL2f(i,j) = sqrt(sum(MEAN_ENERL2) * app.dt * app.PlotIntrvl); %#ok<SAGROW>
            EIGPLOTL2f(:,i,j) = sqrt(sum(EIGPLOTL2,2) * app.dt * app.PlotIntrvl); %#ok<SAGROW>
            YYtL2f(:,i,j) = sqrt(sum(YYtL2 ,2) * app.dt * app.PlotIntrvl); %#ok<SAGROW>

            PLinff(i,j) = max(PLinf); %#ok<SAGROW>
            PiLinff(:,i,j) = max(PiLinf,[],2); %#ok<SAGROW>
            uLinff(i,j) = max(uLinf); %#ok<SAGROW>
            uiLinff(:,i,j) = max(uiLinf,[],2); %#ok<SAGROW>
            vLinff(i,j) = max(vLinf); %#ok<SAGROW>
            viLinff(:,i,j) = max(viLinf,[],2); %#ok<SAGROW>
            rhoLinff(i,j) = max(rhoLinf); %#ok<SAGROW>
            MEAN_ENERLinff(i,j) = max(MEAN_ENERLinf); %#ok<SAGROW>
            EIGPLOTLinff(:,i,j) = max(EIGPLOTLinf,[],2); %#ok<SAGROW>
            YYtLinff(:,i,j) = max(YYtLinf ,[],2); %#ok<SAGROW>
        end
    end
    orderPL2f{j} = log(PL2f(2:end,j) ./ PL2f(1:end-1,j)) / log(2);
    orderPiL2f{j} = (log(PiL2f(:,2:end,j) ./ PiL2f(:,1:end-1,j)) / log(2))';
    orderuL2f{j} = log(uL2f(2:end,j) ./ uL2f(1:end-1,j)) / log(2);
    orderuiL2f{j} = (log(uiL2f(:,2:end,j) ./ uiL2f(:,1:end-1,j)) / log(2))';
    ordervL2f{j} = log(vL2f(2:end,j) ./ vL2f(1:end-1,j)) / log(2);
    orderviL2f{j} = (log(viL2f(:,2:end,j) ./ viL2f(:,1:end-1,j)) / log(2))';
    orderrhoL2f{j} = log(rhoL2f(2:end,j) ./ rhoL2f(1:end-1,j)) / log(2);
    orderMEAN_ENERL2f{j} = log(MEAN_ENERL2f(2:end,j) ./ MEAN_ENERL2f(1:end-1,j)) / log(2);
    orderEIGPLOTL2f{j} = (log(EIGPLOTL2f(:,2:end,j) ./ EIGPLOTL2f(:,1:end-1,j)) / log(2))';
    orderYYtL2f{j} = (log(YYtL2f(:,2:end,j) ./ YYtL2f(:,1:end-1,j)) / log(2))';
    
    orderL2f(:,:,j) = [orderPL2f{j}, orderPiL2f{j}, ...
        orderuL2f{j}, orderuiL2f{j}, ordervL2f{j},orderviL2f{j},...
        orderrhoL2f{j}, orderMEAN_ENERL2f{j}, orderEIGPLOTL2f{j}, orderYYtL2f{j}];
    
    orderPLinff{j} = log(PLinff(2:end,j) ./ PLinff(1:end-1,j)) / log(2);
    orderPiLinff{j} = (log(PiLinff(:,2:end,j) ./ PiLinff(:,1:end-1,j)) / log(2))';
    orderuLinff{j} = log(uLinff(2:end,j) ./ uLinff(1:end-1,j)) / log(2);
    orderuiLinff{j} = (log(uiLinff(:,2:end,j) ./ uiLinff(:,1:end-1,j)) / log(2))';
    ordervLinff{j} = log(vLinff(2:end,j) ./ vLinff(1:end-1,j)) / log(2);
    orderviLinff{j} = (log(viLinff(:,2:end,j) ./ viLinff(:,1:end-1,j)) / log(2))';
    orderrhoLinff{j} = log(rhoLinff(2:end,j) ./ rhoLinff(1:end-1,j)) / log(2);
    orderMEAN_ENERLinff{j} = log(MEAN_ENERLinff(2:end,j) ./ MEAN_ENERLinff(1:end-1,j)) / log(2)';
    orderEIGPLOTLinff{j} = (log(EIGPLOTLinff(:,2:end,j) ./ EIGPLOTLinff(:,1:end-1,j)) / log(2))';
    orderYYtLinff{j} = (log(YYtLinff(:,2:end,j) ./ YYtLinff(:,1:end-1,j)) / log(2))';
    
    orderLinff(:,:,j) = [orderPLinff{j}, orderPiLinff{j}, ...
        orderuLinff{j}, orderuiLinff{j}, ordervLinff{j},orderviLinff{j},...
        orderrhoLinff{j}, orderMEAN_ENERLinff{j}, orderEIGPLOTLinff{j}, orderYYtLinff{j}];
  save(sprintf('%s/Save/%s/RE%d_dt%03d_Final.mat',d0,app.savemat,1/app.nu,1/app.dt))
  
  %Restore mesh resolution for next run 
    app.Nx = app.Nx * 16;
    app.Ny = app.Ny * 16;
end
save(sprintf('%s/Save/%s/RE%d_Final.mat',d0,app.savemat,1/app.nu))
fprintf('Final file save in %s/Save/%s/RE%d_Final.mat',d0,app.savemat,1/app.nu)
