% Lid-Driven Cavity Flow test case setup for DO paper
% Forcing Function
Fu = @(time,p) 0;
Fv = @(time,p) 0;

%Setup boundary size
a    = 1;          b    = 1;
dx   = a./(app.Nx);    dy   = b./(app.Ny); 

%Set boundary conditoins
bcsDP= []; bcsOP= []; bcsNP= [0 0 0 0];        
bcsDu= [0 0 0 1]  ;  bcsNu= [];        
bcsDv= [0 0 0 0]  ;  bcsNv= [];        

%Set number of timesteps
Nt   = ceil(app.T/app.dt);

%% create geometry mask
Nbcs = length(bcsDP) + length(bcsNP);
bnd = [1 1 1 1];
Mask = zeros(app.Ny, app.Nx);

%create Node matrices
[NodeP, Nodeu, Nodev, idsP, idsu, idsv] = NodePad(Mask,bnd,1);

%% Initialize fields
%First grab the coordinates from the reference solution:
d0=fileparts([pwd,filesep]);
d0=d0(1:strfind(d0,'trunk')+4);
if isfield(app,'ref')
    load(sprintf('%s/Save/%s/param.mat',d0,app.ref));
    idsPe = idsP;
    pide = pid;
    XPe = XP;
    YPe = YP;
    XUe =  XU;
    YUe =  YU;
    idsue = idsu;
    uide = uid;
    XVe = XV;
    YVe = YV;
    idsve= idsv;    
end


xp = linspace(dx/2, a-dx/2, app.Nx  );yp = fliplr(linspace(dy/2, b-dy/2, app.Ny  ));
xu = linspace( dx   ,a-dx , app.Nx-1);yu = fliplr(linspace(dy/2, b-dy/2, app.Ny  ));
xv = linspace(dx/2, a-dx/2, app.Nx  );yv = fliplr(linspace(dy    , b-dy, app.Ny-1));
[XP YP] = meshgrid(xp,yp);[XU YU] = meshgrid(xu,yu);[XV YV] = meshgrid(xv,yv);
P = [bcsDP';bcsNP'; 0*idsP];
u = [bcsDu';bcsNu';idsu*0];
v = [bcsDv';bcsNv'; 0*idsv];
ucentid = Nodeu(:,floor(app.Nx/2+1));
vcentid = Nodev(:,floor(app.Nx/2+1));

%%Initialize DO fields:
%Set default number of Monte-Carlo samples
if ~isfield(app,'MC')
    app.MC = 1e5;
end

%Set default variance
if ~isfield(app,'var')
    app.var = 1;
end

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

Nbcs = length(bcsDu) + length(bcsNu);
%%Initialize modes
ui=zeros(length(u),app.S(1));
vi=zeros(length(v),app.S(1));
Pi=zeros(length(P),app.S(1));
var = zeros(app.S(1),1);
% i=1;
% for iy = 1 : ceil(sqrt(app.S(1)))
%     for ix = 1 : iy + 1
%         if (i <= app.S(1))
%             nrm = sqrt((ix^2 + iy^2) * 3) / 8;
%             tmp = - iy * sin(pi * XU * ix).^2 ...
%                 .* cos(pi * YU * iy)  .* sin(pi * YU * iy) ...
%                 / nrm;
%             ui(:,i) = [zeros(Nbcs,1);tmp(idsu)];
%             subplot(121),imagesc(tmp)
%             tmp = ix * sin(pi * YV * iy).^2 ...
%                 .* cos(pi * XV * ix)  .* sin(pi * XV * ix) ...
%                 / nrm;
%             vi(:,i) = [zeros(Nbcs,1);tmp(idsv)];
%             subplot(122),imagesc(tmp)
%             var(i) = app.var * exp(-(ix+iy) + 1);
%             i = i+1;
%             pause
%         end
%     end
% end
% for i = 1 : app.S(1)
%     nrm = sqrt((i^2 + i^2) * 3) / 8;
%     tmp = - i * sin(pi * XU * i).^2 ...
%         .* cos(pi * YU * i)  .* sin(pi * YU * i) ...
%         / nrm;
%     ui(:,i) = [zeros(Nbcs,1);tmp(idsu)];
%     
%     tmp = i * sin(pi * YV * i).^2 ...
%         .* cos(pi * XV * i)  .* sin(pi * XV * i) ...
%         / nrm;
%     vi(:,i) = [zeros(Nbcs,1);tmp(idsv)];
%     
%     var(i) = app.var * exp(-i + 1);    
% end
% uinit = @(x,y)  pi.*(sin(2*pi*y).*(sin(pi*x)).^2);
% vinit = @(x,y) -pi.*(sin(2*pi*x).*(sin(pi*y)).^2);
% for i = 1 : app.S(1)
%     nrm = sqrt(3) * pi / 4;
%     
%     ui(:,i) = [zeros(Nbcs,1);
%         uinit(i * XU(idsu),i * YU(idsu)) / nrm];
%     vi(:,i) = [zeros(Nbcs,1);
%         vinit(i * XV(idsv),i * YV(idsv)) / nrm];
%     
%     var(i) = app.var * exp(-i + 1);    
% end

% uinit = @(x,y)  pi.*(sin(2*pi*y).*(sin(pi*x)).^2);
% vinit = @(x,y) -pi.*(sin(2*pi*x).*(sin(pi*y)).^2);
uinit = @(x,y,M,N) - sin(pi .* x) .* sin(pi .* M .* x)...
    .* (cos(pi .* y) .* sin(pi .* N .* y) + N .* sin(pi .* y) .* cos(pi .* N .* y));
vinit = @(x,y,M,N)   sin(pi .* y) .* sin(pi .* N .* y)...
    .* (cos(pi .* x) .* sin(pi .* M .* x) + M .* sin(pi .* x) .* cos(pi .* M .* x));
i=1;
for iy = 1 : ceil(sqrt(app.S(1)))
    for ix = 1 : iy + 1
        if (i <= app.S(1))
%             nrm = sqrt(3) * pi / 4;
            nrm = sqrt(((ix^2 + 1))/16*(3/2)^(iy==1) +...
                ((iy^2 + 1))/16*(3/2)^(ix==1));
            ui(:,i) = [zeros(Nbcs,1);
                uinit(XU(idsu), YU(idsu), ix, iy) / nrm];
            vi(:,i) = [zeros(Nbcs,1);
                vinit(XV(idsv), YV(idsv), ix, iy) / nrm];
            var(i) = app.var * exp(-(ix+iy) + 1);
            i = i+1;
            pause
            
        end
    end
end

%Initialize stochastic coefficients as normal random with variance as calculated
%in Stoch_Initial_Perturb
YYt=zeros(app.MC,app.S(1));
randn('state', 0);
for i=1:app.S
    YYt(:,i)=normrnd(0,sqrt(var(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
