CylinderSetup_num
app.S = app.S - 1;
ui = sol.ui;
Pi = sol.Pi;
YYt = sol.YYt;
vi = sol.vi;


%% BORROW A LOT FROM THE SOLVER
if ~isfield(app,'IP'),app.IP=[1,1,1];end
if ~exist('bcsDrho','var'),bcsDrho=zeros(size(bcsDP));end
if ~exist('bcsNrho','var'),bcsNrho=zeros(size(bcsNP));end
if ~exist('rho','var'),rho=zeros(size(P));end
if ~exist('rhoi','var'),rhoi=zeros(size(P, 1), app.S);end
if ~isfield(app,'kappa'),app.kappa=0;end
if ~isfield(app,'PrRa')
    app.PrRa=1;
end
if ~exist('bcsOP','var'),bcsOP=[];end
if ~exist('bcsOu','var'),bcsOu=[];end
if ~exist('bcsOv','var'),bcsOv=[];end
if ~isfield(app,'S'),    app.S=0;end
if ~exist('YYt','var'), YYt=0;end
if ~isfield(app,'DOAdapt')
    if length(app.S)==2
        app.DOAdapt=10;
    else
        app.DOAdapt=0;
    end
end
if ~exist('Fcor'), Fcor = @(t,p) zeros(size(p,1),1);end
if ~isfield(app,'DO_cor_iter'), app.DO_cor_iter=1; end


%Find number of boundary conditions
Nbcs = length(bcsDP) + length(bcsNP);

%Check to see how many of the modes have boundary conditions
tmp = (sum(abs(ui(1:Nbcs,:)) + abs(vi(1:Nbcs,:))) > sqrt(eps));
app.Sb = length(find(tmp));
%arrange so that the boundary modes are always first
[tmp, tmp] = sort(tmp, 'descend');
ui = ui(:, tmp);
vi = vi(:, tmp);
Pi = Pi(:, tmp);
if app.Sb > 0
    if ~exist('rhsYibcs', 'var')
        rhsYibcs = @(yi,t) zeros(size(yi));
    end
end

%Set tolerances etc. for GMRES solvers
TOL=sqrt(eps);
MAXIT=20;
RESTARTS=10;

% Low storage Runge-Kutta coefficients
rk4a = [            0.0 ...
    -567301805773.0/1357537059087.0 ...
    -2404267990393.0/2016746695238.0 ...
    -3550918686646.0/2091501179385.0  ...
    -1275806237668.0/842570457699.0];
rk4b = [ 1432997174477.0/9575080441755.0 ...
    5161836677717.0/13612068292357.0 ...
    1720146321549.0/2090206949498.0  ...
    3134564353537.0/4481467310338.0  ...
    2277821191437.0/14882151754819.0];
rk4c = [             0.0  ...
    1432997174477.0/9575080441755.0 ...
    2526269341429.0/6820363962896.0 ...
    2006345519317.0/3224310063776.0 ...
    2802321613138.0/2924317926251.0];

%% Constructing A and RHS
[Arho rhsrho] =  mk_DiffusionOperator(app.Nx, app.Ny, dx, dy, bcsDrho,bcsNrho, NodeP,'P');
[AP rhsPbc] =  mk_DiffusionOperator(app.Nx, app.Ny, dx, dy, bcsDP,bcsNP, NodeP,'P');
[Au rhsubc] =  mk_DiffusionOperator(app.Nx, app.Ny, dx, dy, bcsDu,bcsNu, Nodeu,'u');
[Av rhsvbc] =  mk_DiffusionOperator(app.Nx, app.Ny, dx, dy, bcsDv,bcsNv, Nodev,'v');
[AP rhsPbc] = FixDiffPObcs(AP,bcsOP,Nbcs,...
    rhsPbc,NodeP,dx,dy);bcsDP(bcsOP)=[];
[Au rhsubc] = FixDiffPObcs(Au,bcsOu,Nbcs,...
    rhsubc,Nodeu,dx,dy);bcsDu(bcsOu)=[];
[Av rhsvbc] = FixDiffPObcs(Av,[],Nbcs,...
    rhsvbc,Nodev,dx,dy);bcsDv([])=[];
[Au Av rhsubc rhsvbc] = FixDiffuvDbcs(Au,Av,bcsDu,bcsDv,...
    rhsubc,rhsvbc,Nodeu, Nodev, dx,dy,Nbcs);

if (sum(abs(Fcor(0, [XV(idsv) YV(idsv)])))>0)
    [Au2v Av2u] = CoriolisAVG(app.Nx, app.Ny, Nodeu, Nodev, idsu, idsv, Nbcs);
    Au2v = spdiags(Fcor(0 + app.dt, [XV(idsv) YV(idsv)]), 0, ...
        length(idsv), length(idsv)) * Au2v;
    Av2u = spdiags(Fcor(0 + app.dt, [XU(idsu) YU(idsu)]), 0, ...
        length(idsu), length(idsu)) * Av2u;
    %Au2v = sparse(diag(Fcor(0 + app.dt,[XV(idsv) YV(idsv)])))*Au2v;
    %Av2u = sparse(diag(Fcor(0 + app.dt,[XU(idsu) YU(idsu)])))*Av2u;
else
    Au2v = spalloc(size(v,1) - Nbcs,size(u,1),0);
    Av2u = spalloc(size(u,1) - Nbcs,size(v,1),0);
end

AU=-app.nu*Au;
AV=-app.nu*Av;
ARHO = -app.kappa*Arho;

Au = app.nu*Au + 1/app.dt*speye(size(Au));
Av = app.nu*Av + 1/app.dt*speye(size(Av));
Arho = app.kappa*Arho + 1/app.dt*speye(size(Arho));

%% LU factorization of matrices
imposemean = false;
if isempty(bcsDP)
    imposemean = true;
end
%First order matrix to give best spartsity pattern
rhoperm = symamd(Arho);
Pperm = symamd(AP);
uperm = symamd(Au);
vperm = symamd(Av);

AP=AP(Pperm,Pperm);

% % LU factorization
[Lrho Urho] = lu(Arho(rhoperm,rhoperm));
[LP UP] = lu(AP);
[Lu Uu] = lu(Au(uperm,uperm));
[Lv Uv] = lu(Av(vperm,vperm));

% Get interior ids
pid = Nbcs+1:length(P);
uid = Nbcs+1:length(u);
vid = Nbcs+1:length(v);
srcid     = 2:app.Nx+1;
srcTopid  = 2:app.Ny;
srcBotid  = 3:app.Ny+1;
Neu = [length(bcsDu)+1,Nbcs];
Nev = [length(bcsDv)+1,Nbcs];
NbcsDrho = length(bcsDrho); %Needed for advection
NbcsDu = length(bcsDu); %Needed for advection
NbcsDv = length(bcsDv); %Needed for advection

%% Initialize variables and sizes
frho = zeros(length(pid),1);
fu = zeros(length(uid),1);
fv = zeros(length(vid),1);

if app.S(1)>0 || length(app.S)>1
    frhoi = zeros(length(pid),app.S(1));
    fui = zeros(length(uid),app.S(1));
    fvi = zeros(length(vid),app.S(1));
    frhoij = zeros(length(pid),app.S(1)^2);
    fuij= zeros(length(uid),app.S(1)^2);
    fvij= zeros(length(vid),app.S(1)^2);
else
    frhoi = 0;
    frhoij = 0;
    fui = 0;
    fvi = 0;
    fuij= 0;
    fvij= 0;
end

%If stochastic modes were not initialized in Setup script,
%initialize automatically.
DOInit;
if ~exist('rhoi','var') && (app.S(1)>0 || length(app.S)>1)
    rhoi = zeros(length(rho),app.S(1));
end

%% Solving system
q = zeros(size(P));
CYY=cov(YYt);

T=app.dt*sol.k;
%Need to save ui00 etc to do the projections correctly
ui00 = ui;
vi00 = vi;
rhoi00 = rhoi;

%Calculate statistics
MY = mean(YYt);
for i=1:app.S(1)
    YYt(:,i) = YYt(:,i) - MY(i);
end
CYY=cov(YYt);
CYYINV = pinv(CYY, max(CYY(:)) * sqrt(eps));
MYYY=mom_thrd(YYt);
if app.S(1)>0
    MYYY=reshape(MYYY,app.S(1)*app.S(1),app.S(1));
end
if app.S(1) ==0
    CYY=0;
end
%Solve mean field
[dpdx dpdy]   = Grad2D(P,NodeP,dx,dy,idsu,idsv);
% Calculate RHS vectors for mean u, v, and rho
[fu fv] = UVadvect_TVD_mex(app.Nx, app.Ny, dx, dy, app.dt, u, v, u, v,Nodeu, Nodev, idsu, idsv, Nbcs, NbcsDu, NbcsDv);
fu = fu + Av2u*v;
src = app.PrRa * 0.5*( rho(NodeP(srcTopid,srcid)) + rho(NodeP(srcBotid,srcid)) );
src = src(idsv);
fv = fv - Au2v*u -src;
frho = SCAadvect_TVD_mex(app.Nx, app.Ny, dx, dy, app.dt, rho, u, v, NodeP, Nodeu, Nodev, idsP,Nbcs,NbcsDrho,NbcsDu,NbcsDv);

%calculate fui and fvi forcings, as well as fuij, and fvij
for i=1:app.S(1)
    [fui(:,i) fvi(:,i)] = UVadvect_TVD_mex(app.Nx, app.Ny, dx, dy, app.dt,...
        ui(:,i), vi(:,i),u      , v      , ...
        Nodeu, Nodev, idsu, idsv, Nbcs, NbcsDu, NbcsDv);
    [tmp tmp1]          = UVadvect_TVD_mex(app.Nx, app.Ny, dx, dy, app.dt,...
        u      , v      ,ui(:,i), vi(:,i),...
        Nodeu, Nodev, idsu, idsv, Nbcs, NbcsDu, NbcsDv);
    [tmp2 tmp3]          = UVadvect_TVD_mex(app.Nx, app.Ny, dx, dy, app.dt,...
        u      , v      ,-ui(:,i), -vi(:,i),...
        Nodeu, Nodev, idsu, idsv, Nbcs, NbcsDu, NbcsDv);
    
    src = app.PrRa * 0.5*( rhoi(NodeP(srcTopid,srcid),i) + rhoi(NodeP(srcBotid,srcid),i));
    src = src(idsv);
    fui(:,i) = fui(:,i) + 0.5*(tmp - tmp2)  + Av2u*vi(:,i);
    fvi(:,i) = fvi(:,i) + 0.5*(tmp1 - tmp3) - Au2v*ui(:,i) - src;
    
    frhoi(:,i) = SCAadvect_TVD_mex(app.Nx, app.Ny, dx, dy, app.dt, rhoi(:,i),...
        u, v, NodeP, Nodeu, Nodev, idsP,Nbcs,NbcsDrho,NbcsDu,NbcsDv) ...
        + 0.5 * SCAadvect_TVD_mex(app.Nx, app.Ny, dx, dy, app.dt, rho,...
        ui(:,i), vi(:,i), NodeP, Nodeu, Nodev, idsP,Nbcs,NbcsDrho,NbcsDu,NbcsDv) ...
        - 0.5 * SCAadvect_TVD_mex(app.Nx, app.Ny, dx, dy, app.dt, rho,...
        -ui(:,i), -vi(:,i), NodeP, Nodeu, Nodev, idsP,Nbcs,NbcsDrho,NbcsDu,NbcsDv);
    %For stochastic boundary forcing, need to change above two lines
    for j=1:app.S(1)
        [fuij(:,i+(j-1)*app.S(1)) fvij(:,i+(j-1)*app.S(1))] = ...
            UVadvect_TVD_mex(app.Nx, app.Ny, dx, dy, app.dt,...
            ui(:,i),vi(:,i),ui(:,j),vi(:,j), ...
            Nodeu, Nodev, idsu, idsv, Nbcs, NbcsDu, NbcsDv);
        [tmp tmp1] = ...
            UVadvect_TVD_mex(app.Nx, app.Ny, dx, dy, app.dt,...
            ui(:,i),vi(:,i),-ui(:,j),-vi(:,j), ...
            Nodeu, Nodev, idsu, idsv, Nbcs, NbcsDu, NbcsDv);
        frhoij(:,i+(j-1)*app.S(1)) = 0.5 * SCAadvect_TVD_mex(app.Nx, app.Ny, dx, dy, app.dt, rhoi(:,i),...
            ui(:,j), vi(:,j), NodeP, Nodeu, Nodev, idsP,Nbcs,NbcsDrho,NbcsDu,NbcsDv)...
            - 0.5 * SCAadvect_TVD_mex(app.Nx, app.Ny, dx, dy, app.dt, rhoi(:,i),...
            -ui(:,j), -vi(:,j), NodeP, Nodeu, Nodev, idsP,Nbcs,NbcsDrho,NbcsDu,NbcsDv);
        %%To be consistent with the notation in the paper, let us take
        %%the negative sign of the above expression:
        frhoij(:,i+(j-1)*app.S(1)) = -frhoij(:,i+(j-1)*app.S(1));
        fuij(:,i+(j-1)*app.S(1)) = 0.5*(tmp - fuij(:,i+(j-1)*app.S(1)));
        fvij(:,i+(j-1)*app.S(1)) =  0.5*(tmp1 - fvij(:,i+(j-1)*app.S(1)));
    end
end


%Advance Yi's using RK4
if app.S(1)>0
    fui_diffusion = fui+AU*ui(uid,:);
    fvi_diffusion = fvi+AV*vi(vid,:);
    frhoi_diffusion = frhoi + ARHO*rhoi(pid,:);
end
RHSY = zeros(size(YYt));
resYYt = zeros(size(YYt));
Tlocal0 = T + app.dt * rk4c(1);
%Calculate YiYj's
for i=1:app.S(1)
    for j=1:app.S(1)
        YiYj(:, i + (j - 1) * app.S(1)) = ...
            YYt(: , i) .* YYt(:, j);
    end
end

RHSY(:, 1:app.Sb) = rhsYibcs(YYt(:, 1:app.Sb), Tlocal0);
for i=app.Sb+1:app.S(1)
    A1j =  dx * dy * (app.IP(1) * (fui_diffusion)' * ui(uid, i) + ...
        app.IP(2) * (fvi_diffusion)' * vi(vid, i) + ...
        app.IP(3) * (frhoi_diffusion)' * rhoi(pid, i)); %column vector
    A2j =  dx * dy * (app.IP(1) * fuij' * ui(uid, i) + ...
        app.IP(2) * fvij' * vi(vid, i) + ...
        app.IP(3) * frhoij' * rhoi(pid, i));
    RHSY(:, i) = YYt * A1j + CYY(:)' * A2j(:) - YiYj * A2j(:);
end