% Lid-Driven Cavity Flow test case

%First call the setup file from the LidDrivenCavityFlow, and add the DO part.
LidDrivenCavitySetup;

%%Initialize DO fields:
%Because it takes long to do the initialization, we save the results for the 
%next time.
d0=fileparts([pwd,filesep]);
d0=d0(1:strfind(d0,'trunk')+4);

if exist(sprintf('%s/Save/DO_ICs/uivi_%d_%d.mat',d0,app.Nx,app.Ny),'file')
    eval(sprintf('load %s/Save/DO_ICs/uivi_%d_%d.mat',d0,app.Nx,app.Ny));
    if isfield(app,'var')
        var = var/var(1)*app.var(1);
    else
        var = var/var(1)*1e-6;
    end
else
    [ui, vi, var] = Stoch_Initial_Perturb(XU,YU,XV,YV,Nbcs,...
        length(bcsDP)+length(bcsOP),NodeP,12);
    if isfield(app,'var')
        var = var*app.var(1);
    else
        var = var*1e-6;
    end
    display('Created initial conditions for DO simulation.')
    display('Paused, press any key to continue. Will save upon continue.')
    pause;
    eval(sprintf('save %s/Save/DO_ICs/uivi_%d_%d.mat ui vi var',d0,app.Nx,app.Ny));
    close all;
end

%Set default number of Monte-Carlo samples
if ~isfield(app,'MC')
    app.MC=4000;
end

app.S=min(app.S(1),12);

%%Initialize modes
if app.S(1) < 6
    %For small S, these modes seemed to be a better choice for initialization
    ui=ui(:,[1,4,5,6,7]);
    vi=vi(:,[1,4,5,6,7]);
    I=[1,4,5,6,7];
end
ui=ui(:,1:1:app.S(1));
vi=vi(:,1:1:app.S(1));
Pi=zeros(length(P),app.S(1));

%Initialize stochastic coefficients as normal random with variance as calculated
%in Stoch_Initial_Perturb
YYt=zeros(app.MC,app.S(1));
for i=1:app.S
    YYt(:,i)=normrnd(0,sqrt(var(I(i))),app.MC,1);
end

%Need these for Orthonormalization
uid = Nbcs+1:length(u);
vid = Nbcs+1:length(v);
%Ortho-Normalization of basis fields to avoid round-off errors
for i=1:app.S(1)
    %Commented below in order to check with Themis' code
     for j=1:app.S(1)
         if i~=j
         coef = InProd(ui(:,i),ui(:,j),dx,dy,uid) + ...
                InProd(vi(:,i),vi(:,j),dx,dy,vid);
         ui(:,i) = ui(:,i) - coef * ui(:,j);
         vi(:,i) = vi(:,i) - coef * vi(:,j);
         end
     end
    nrm=sqrt(InProd(ui(:,i),ui(:,i),dx,dy,uid)+InProd(vi(:,i),vi(:,i),dx,dy,vid));
    ui(:,i) = ui(:,i)/nrm;
    vi(:,i) = vi(:,i)/nrm;
end
