% Naviers Stokes Flow Lock Exchange test case

LockExchangeSetup;

%Set default number of Monte-Carlo samples
if ~isfield(app, 'MC'), error('app.MC not specified'); end;
if ~isfield(app, 'S'), error('app.S not specified'); end;
if app.S(1)>2
    display('Only works for S=1,2');
    app.S(1)=2;
end
%%Initialize modes
ui=zeros(length(u),app.S(1));
vi=zeros(length(v),app.S(1));
Pi=zeros(length(P),app.S(1));
rhoi=ones(length(rho),app.S(1));
if isfield(app,'var')
    var = exp(0:-1:-app.S(1)+1)*app.var;
else
    var = exp(0:-1:-app.S(1)+1)*1e-16;
end

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

%Need these for Orthonormalization
uid = Nbcs+1:length(u);
vid = Nbcs+1:length(v);
pid = Nbcs+1:length(rho);

% Now we add a stochastic initial rho!
if isfield(app,'Stochvar')
    var(1) = app.Stochvar;
else
    var(1) = 1e-16;
end

%Now initialize density
if app.S(1)>0
    rhoi(pid,1) = -rho(pid);
    Yi(:,1)=9.81/2*[ones(app.MC/2,1);-ones(app.MC/2,1) ];
    if app.S(1)>1
         Yi(:,1)=[9.81/4*ones(app.MC/4,1);9.81/4*3*ones(app.MC/4,1);
                  -9.81/2*ones(app.MC/2,1)];     
    end
    Yi = Yi * 0.5;
else
    app.S(1)=1;
    Yi=zeros(app.MC,1);
    ui=zeros(length(u),app.S(1));
    vi=zeros(length(v),app.S(1));
    Pi=zeros(length(P),app.S(1));
    rhoi=ones(length(rho),app.S(1));
end
for i=3:app.S(1)
    tmp = sin(XP*pi).*sin(YP*pi);
    rhoi(pid,i)=tmp(idsP);
end

%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) + ...
                InProd(rhoi(:,i),rhoi(:,j),dx,dy,pid);
            ui(:,i) = ui(:,i) - coef * ui(:,j);
            vi(:,i) = vi(:,i) - coef * vi(:,j);
            rhoi(:,i) = rhoi(:,i) - coef * rhoi(:,j);
        end
    end
    nrm=sqrt(InProd(ui(:,i),ui(:,i),dx,dy,uid) + ...
        InProd(vi(:,i),vi(:,i),dx,dy,vid) + ...
        InProd(rhoi(:,i),rhoi(:,i),dx,dy,pid));
    ui(:,i) = ui(:,i)/nrm;
    vi(:,i) = vi(:,i)/nrm;
    rhoi(:,i) = rhoi(:,i)/nrm;
end
