%clear all, clc, clf, close all
function [] = LockExchangeRunScript(Gr,runref)
if nargin < 1 || isempty(Gr)
    Gr = 1.25e6;
end
if nargin < 2 || isempty(runref)
    runref = true;
end

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

%Create the app for the Reference solution
sscriptref = 'LockExchangeSetup_numdet';
pscriptref = 'GenericSaveScriptRho';
app.Nx = 256;
app.Ny = 256;
app.dt = 1 / app.Nx / 4 / 2;
app.T = 5;
app.PlotIntrvl = 64 * 2;
app.Gr = Gr;
app.Sc = 1;

%app.saveimg = app.savemat;
app.var = 1 / 64;
drho = [0.5, 0.42, 0.37, 0.31] * 2;
% drho = [0.5, 0.42, 0.38, 0.3] * 2;
% drho = [0.5, 0.42, 0.38] * 2;
% drho = [0.5,0.3]*2;

%Run the reference solution
if runref
    for i = 1:length(drho)
        app.drho = drho(i);
        app.savemat = ['DOnumPap/LockRefDet_',num2str(Gr),'_',...
            num2str(app.Nx),'_',num2str(app.drho*100)];
        NavierStokesInc(app,sscriptref,pscriptref);
%         NavierStokesInc(app,sscript,pscript);
        %and noisy deterministic runs
        app.noise = 1e-1;
        app.drho = drho(i) + 1e-1;
        app.savemat = ['DOnumPap/LockRefDetNoise_',num2str(Gr),'_',...
            num2str(app.Nx),'_',num2str(app.drho*100)];
        NavierStokesInc(app,sscriptref,pscriptref);
        app = rmfield(app,'noise');
    end
end
app.S = 3;
app.MC = 4;
app.drho = mean(drho);

%Run the stochastic case
if isfield(app,'saveimg')
    app = rmfield(app,'saveimg');
end
app.ref = ['DOnumPap/LockRefDet_',num2str(Gr),'_',num2str(app.Nx),'_'];
app.ref_noise = ['DOnumPap/LockRefDetNoise_',num2str(Gr),'_',num2str(app.Nx),'_'];
app.dref = drho * 100;
app.savemat = ['DOnumPap/LockStoch_',num2str(Gr),'_NX',num2str(app.Nx)]
% app.PlotIntrvl = 1
% app.saveimg = app.savemat;
% app.dt = app.dt / 2;
% app.PlotIntrvl = app.PlotIntrvl * 2;
%This uses the unbiased estimator for the CYY, which is wrong
% NavierStokesInc_DO_rho(app,sscript,pscript);
%This uses the biased estimator:
app.torder = 1;
NavierStokesInc_DO_rho_DIRK_cov2_orthofix(app,sscript,pscript);
% NavierStokesInc_DO_rho_DIRK_cov2(app,sscript,pscript);

app.savemat = ['DOnumPap/LockStoch_GSorth_',num2str(Gr),'_NX',num2str(app.Nx)]
app.GSorth = []
%This uses the unbiased estimator for the CYY, which is wrong
% NavierStokesInc_DO_rho(app,sscript,pscript);
%This uses the biased estimator:
app.torder = 1;
NavierStokesInc_DO_rho_DIRK_cov2_orthofix(app,sscript,pscript);
%NavierStokesInc_DO_rho_DIRK_cov2(app,sscript,pscript);


% %higher order DO attempt:
% % app.torder = 2;
% % NavierStokesInc_DO_rho_DIRK_cov2(app,sscript,pscript);
% % NavierStokesInc_DO_rho_tm2(app,sscript,pscript);

% %% Run only 2 cases
% app.S = 1
% app.dref = app.dref([1,4]);
% app.drho = mean(drho([1 4]));
% NavierStokesInc_DO_rho(app,sscript,pscript);
% % NavierStokesInc_DO_rho1(app,sscript,pscript);

% app.k = app.PlotIntrvl * 8;
% sscript = 'LockExchangeSetup_num3';
% NavierStokesInc_DO_rho(app,sscript,pscript);
