% Stochastic Rayleigh-Benard test case
% (Navier-Stokes with Boussinesq Approximation)

RayleighBenardSetup;

% User specifies the functions that provide initial values at given coordinates.

a1=2.365; a2=3.9266;
P_init   = @(X, Y) zeros(size(X));
u_init   = @(X, Y) ...
    0.3 * (cos(a1*(X-1/2))./cos(a1/2) - cosh(a1*(X-0.5))./cosh(a1/2)) ...
    .* (-a1*sin(a1*(Y-1/2))./cos(a1/2) - a1*sinh(a1*(Y-0.5))./cosh(a1/2)) ...
    + 0.2 * (sin(a2*(X-1/2))./sin(a2/2) - sinh(a2*(X-0.5))./sinh(a2/2)) ...
    .* (a2*cos(a2*(Y-1/2))./sin(a2/2) - a2*cosh(a2*(Y-0.5))./sinh(a2/2));
v_init   = @(X, Y) ...
    -0.3 * (cos(a1*(Y-1/2))./cos(a1/2) - cosh(a1*(Y-0.5))./cosh(a1/2)) ...
    .* (-a1*sin(a1*(X-1/2))./cos(a1/2) - a1*sinh(a1*(X-0.5))./cosh(a1/2)) ...
    - 0.2 * (sin(a2*(Y-1/2))./sin(a2/2) - sinh(a2*(Y-0.5))./sinh(a2/2)) ...
    .* (a2*cos(a2*(X-1/2))./sin(a2/2) - a2*cosh(a2*(X-0.5))./sinh(a2/2));
rho_init = @(X, Y) 1-Y;


P   = [bcsDq; bcsNq; P_init(XP(idsP), YP(idsP))];
u   = [bcsDu; bcsNu; u_init(XU(idsu), YU(idsu))];
v   = [bcsDv; bcsNv; v_init(XV(idsv), YV(idsv))];
rho = [bcsDrho; bcsNrho; rho_init(XP(idsP), YP(idsP))];



%% modes definition

ui=zeros(length(u),app.S(1));
vi=zeros(length(v),app.S(1));
Pi=zeros(length(P),app.S(1));
rhoi=zeros(length(rho),app.S(1));

for i=1:app.S(1)
    ui(:,i)=u;
    vi(:,i)=v;
end

for i = 1:app.S(1)
    rhomod1=@(X,Y) 2*cos(i*pi*X).*sin(pi*Y);
    rhoi(Nbcs+1:length(rho),i)=rhomod1(XP(idsP),YP(idsP))';
end

app.MC=5*10^4;

Yi=zeros(app.MC,app.S(1));

vv=[2*10^-1, 1*10^-1, 3*10^-7, 4*10^-7, 5*10^-7]; 
for i = 1:app.S(1)
    zz=(normrnd(0,sqrt(vv(i)),app.MC/2,1));
    Yi(:,i) = [zz;-zz]'; 
end

% Get interior ids
pid = Nbcs+1:length(P);
uid = Nbcs+1:length(u);
vid = Nbcs+1:length(v);

[Yi, ui, vi, Pi, rhoi]...
        = DOorthnorm(app, dx, dy, Yi, ui, uid, vi, vid, Pi, pid, rhoi);