function [P u v rho NodeP Nodeu Nodev] =...
    NavierStokesInc_DO_mx(app,SetupScript,PlotScript)
% function [P u v rho NodeP Nodeu Nodev] =
% NavierStokesInc_DO(app,SetupScript,PlotScript)
% Solves the Boussinesq equations using a simple Finite Volume Incremental
% Pressure Correction in Rotational form projection method.
% This function uses a first-order implicit Euler integration scheme, and
% uses LU factorization to solve the matrices.The advection uses a TVD
% scheme when the mean is involved, and an averaged TVD scheme when
% advected with the modes invovled. Advection routines are mexed.
%
%INPUTS:
%       app:	    The application data structured. At a minimum it should
%		     contain the following fields (which may be added through
%                    the SetupScript).
%          app.Nx:         Number of interior nodes in x-direction (integer number)
%          app.Ny:         Number of interior nodes in y-direction (integer number)
%          app.T:          Final time
%          app.dt:         Timestep size
%          app.nu:         Viscosity
%          app.kappa:      Diffusivity of density (optional)
%          app.S:          Size of the stochastic dimension. If length(app.S)==1.
%                           If length(app.S)==2, then we have adaptive number of
%                           stochastic dimensions:
%             app.S(1):       Minimum number of stochastic dimensions
%             app.S(2):       Maximum number of stochastic dimensions
%          app.var           Variance of largest mode at initialization
%                             (Optional) for DO solvers
%          app.addvar        If smallest variance mode has this variance,
%                              a new mode is added for adaptive DO solvers
%                             (Optional)
%          app.subbvar       If smallest variance mode has this variance, and
%                              second smallest variance is less than twice
%                              this variance, then the smallest variance mode
%                              is subtraced for adaptive DO solvers (Optional)
%          app.MC:         Number of monte-carlo samples for stochastic
%                           coefficients (Default = 10,000)
%          app.DOAdapt:    Number of time integration intervals between checking
%                           if the number of stochastic modes should be adjusted
%                           for adaptive stochastic dimension case
%                           (length(app.S)==2) Default value is 10.
%          app.PlotIntrvl: Number of time integration intervals between calling
%                           "plotscript."
%       PlotScript: (String input) of filename used to plot specific problem.
%                   Can also be used to save outputs and perform other tasks.
%                       Note: Variables u, v, and P, will be available for
%                       plotting.
%       SetupScript: String input of filename used to setup the problem.
%           This file needs to correctly create the following variables:
%	    Nt:     Number of timesteps (integer number)
%           dx:     Size of control volume in x-direction
%           dy:     Size of control volume in y-direction
%           bcsDrho:Density Dirichlet boundary conditions (optional)
%           bcsDP:  Pressure Dirichlet boundary conditions
%           bcsDu:  U-velocity Dirichlet boundary conditions
%           bcsDv:  V-velocity Dirichlet boundary conditions
%           bcsNrho:Density Neumann boundary conditions (optional)
%           bcsNP:  Pressure Neumann boundary conditions
%           bcsNu:  U-velocity Neumann boundary conditions
%           bcsNv:  V-velocity Neumann boundary conditions
%           NodeP:  Pressure node-numbering matrix (see NodePad.m)
%           Nodeu:  U-velocity node-numbering matrix (see NodePad.m)
%           Nodev:  V-velocity node-numbering matrix (see NodePad.m)
%           idsP:   Interior Pressure node id's (see NodePad.m)
%           idsu:   Interior U-velocity node id's (see NodePad.m)
%           idsv:   Interior V-velocity node id's (see NodePad.m)
%           Fu(time,p):   Forcing function for U-velocity
%                         p =[ XU(idsu),YU(idsu) ]
%           Fv(time,p):   Forcing function for V-velocity
%                         p =[ XV(idsv),YV(idsv) ]
%           XU:     Matrix of x-coordinates corresponding to U-velocity CV
%                    centers. XU(idsu) used as 'x' input to Fu. (see meshgrid)
%           XV:     Matrix of x-coordinates corresponding to V-velocity CV
%                    centers. YU(idsu) used as 'y' input to Fv. (see meshgrid)
%           YU:     Matrix of y-coordinates corresponding to U-velocity CV
%                    centers. XV(idsv) used as 'x' input to Fu. (see meshgrid)
%           YV:     Matrix of y-coordinates corresponding to V-velocity CV
%                    centers. YV(idsv) used as 'y' input to Fv. (see meshgrid)
%           P:      Initial Pressure vector. The first #Nbc terms contain
%                   the boundary condition values.
%           u:      Initial U-velocity vector. The first #Nbc terms contain
%                   the boundary condition values.
%           v:      Initial V-velocity vector. The first #Nbc terms contain
%                   the boundary condition values.
%
%OUTPUTS:
%        P:      Final Pressure vector, including the boundary
%                 conditions.
%        u:      Final U-velocity vector, including the boundary
%                 conditions.
%        v:      Final V-velocity vector, including the boundary
%                 conditions.
%        rho:    Final density vector, including the boundary
%                 conditions.
%        NodeP:  Pressure node-numbering matrix (useful for plotting)
%        Nodeu:  U-velocity node-numbering matrix (useful for plotting)
%        Nodev:  V-velocity node-numbering matrix (useful for plotting)
%
%   Uses Matlab functions: symamd, lu
%   See also: mk_DiffusionOperator, FixDiffuvDbcs, Grad2D, Div2D
%
% Authors:   Matt Ueckermann, Themis Sapsis,
%            Pierre Lermusiaux, Pat Haley for MIT course 2.29

%Run setup script
eval(SetupScript)
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 ~isfield(app,'kappa'),app.kappa=0;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 ~isfield(app,'DO_cor_iter'), app.DO_cor_iter=1; end

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

%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,bcsOv,Nbcs,...
    rhsvbc,Nodev,dx,dy);bcsDv(bcsOv)=[];
[Au Av rhsubc rhsvbc] = FixDiffuvDbcs(Au,Av,bcsDu,bcsDv,...
    rhsubc,rhsvbc,Nodeu, Nodev, dx,dy,Nbcs);
Arho = app.kappa*Arho + 1/app.dt*speye(size(Arho));

AU=-app.nu*Au;
AV=-app.nu*Av;

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

%% 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];
NbcsDP = length(bcsDrho); %Needed for advection
NbcsDu = length(bcsDu); %Needed for advection
NbcsDv = length(bcsDv); %Needed for advection

%% Initialize variables and sizes
fu = zeros(length(uid),1);
fv = zeros(length(vid),1);
if app.S(1)>0 || length(app.S)>1
    fui = zeros(length(uid),app.S(1));
    fvi = zeros(length(vid),app.S(1));
    fuij= zeros(length(uid),app.S(1)^2);
    fvij= zeros(length(vid),app.S(1)^2);
else
    fui = 0;
    fvi = 0;
    fuij= 0;
    fvij= 0;
end

%If stochastic modes were not initialized in Setup script,
%initialize automatically.
DOInit;

%% Solving system
q = zeros(size(P));
CYY=cov(YYt);
T=0;k=0;  eval(PlotScript);
if app.DOAdapt > 0
    DOAdapt_mex;
end
display('Paused at initial condition. Hit any key to continue');drawnow;pause
for k=1:Nt
    T=app.dt*k;
    
    %Need to save ui00 etc to do the projections correctly
    ui00 = ui;
    vi00 = vi;    
    
    %Calculate statistics
    MY=mean(YYt);
    for i=1:app.S
        YYt(:,i) = YYt(:,i) - MY(i);
        %         u(uid) = u(uid) + MY(i) * ui(uid,i);
        %         v(vid) = v(vid) + MY(i) * vi(vid,i);
        %         P(pid) = P(pid) + MY(i) * Pi(pid,i);
    end
    CYY=cov(YYt);
    CYYINV = inv(CYY);
    MYYY=mom_thrd(YYt);
    if app.S>0
        MYYY=reshape(MYYY,app.S*app.S,app.S);
    end
    
    %Solve mean field
    [dpdx dpdy]   = Grad2D(P,NodeP,dx,dy,idsu,idsv);
    % Calculate RHS vectors for u and v
    [fu fv] = UVadvect_TVD_mex(app.Nx, app.Ny, dx, dy, app.dt, u, v, u, v,Nodeu, Nodev, idsu, idsv, Nbcs, NbcsDu, NbcsDv);
    
    %calculate fui and fvi forcings, as well as fuij, and fvij
    for i=1:app.S
        [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);
        
        fui(:,i) = fui(:,i) + 0.5*(tmp - tmp2);
        fvi(:,i) = fvi(:,i) + 0.5*(tmp1 - tmp3);
        
        
        %For stochastic boundary forcing, need to change above two lines
        for j=1:app.S
            [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);
            %%To be consistent with the notation in the paper, let us take
            %%the negative sign of the above expression:
            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>0
        fui_diffusion = fui+AU*ui(uid,:);
        fvi_diffusion = fvi+AV*vi(vid,:);
    end
        RHSY=zeros(size(YYt));
    %Calculate YiYj's
     YYt0 = YYt;
     for irk = 1:2
         %RK time stepping
         YYt = YYt0 + app.dt * RHSY;
         for i=1:app.S
             for j=1:app.S
                 YiYj(:,i+(j-1)*app.S)=YYt(:,i).*YYt(:,j);
             end
         end
         for i=1:app.S
             A1j =  dx*dy*((fui_diffusion)'*ui(uid,i)+...
                 (fvi_diffusion)'*vi(vid,i)); %column vector
             A2j =  dx*dy*(fuij'*ui(uid,i)+fvij'*vi(vid,i));
             RHSY(:,i) = YYt*A1j+CYY(:)'*A2j(:)-YiYj*A2j(:);
         end
     end
    YYt = YYt0 + 0.5 * (app.dt * RHSY + YYt - YYt0);
   
    
    % Advance mean field
    src = 0.5*( rho(NodeP(srcTopid,srcid)) + rho(NodeP(srcBotid,srcid)) );
    src = src(idsv);
    RHSu = 1/app.dt*u(uid) + fu - fuij*CYY(:) + rhsubc*app.nu - dpdx + ...
        Fu(T + app.dt,[XU(idsu) YU(idsu)]);
    RHSv = 1/app.dt*v(vid) + fv - fvij*CYY(:) + rhsvbc*app.nu - dpdy +...
        Fv(T + app.dt,[XV(idsv) YV(idsv)])-src;
    
    u(uid(uperm))=Uu\(Lu\RHSu(uperm));
    v(vid(vperm))=Uv\(Lv\RHSv(vperm));
    
    fuijmyyy = fuij*MYYY;
    fvijmyyy = fvij*MYYY;
    %Advance ui's and vi's
    for i=1:app.S
        f1u = fui(:,i) - fuijmyyy*CYYINV(:,i);
        f1v = fvi(:,i) - fvijmyyy*CYYINV(:,i);
        f2u = f1u + AU*ui(uid,i);
        f2v = f1v + AV*vi(vid,i);
        %Calculate inner product terms <ffiUj> -- there has to be app.S of them
        InProdi=zeros(app.S,1);
        for j=1:app.S
            InProdi(j) = (f2u'*ui00(uid,j)+f2v'*vi00(vid,j))*dx*dy;
        end
        
        %We're going to treat the Laplacian for THIS mode implicitly.
        %but the inner product treats the Laplcian EXplicitly
        [dpdx dpdy]   = Grad2D(Pi(:,i),NodeP,dx,dy,idsu,idsv);
        RHSu = 1/app.dt*ui(uid,i) + f1u - dpdx...
            -ui00(uid,:)*InProdi;
        %NOTE! I am assuming that ui's are zero on the boundary. To
        %introduce stochastic boundary forcings, add to the above line.
        RHSv = 1/app.dt*vi(vid,i) + f1v - dpdy...
            -vi00(vid,:)*InProdi;
        %NOTE! I am assuming that vi's are zero on the boundary. To
        %introduce stochastic boundary forcings, add to the above line.
        ui(uid(uperm),i)=Uu\(Lu\RHSu(uperm));
        vi(vid(vperm),i)=Uv\(Lv\RHSv(vperm));
    end
    
    %Calculate RHS vector for P
    divU = Div2D(u,v,Nodeu,Nodev,dx,dy,idsP,Neu,Nev);
    RHSq = rhsPbc - 1/app.dt*divU;
    q(pid(Pperm))=UP\(LP\RHSq(Pperm));
    
    %Correct Velocity
    [dqdx dqdy]   = Grad2D(q,NodeP,dx,dy,idsu,idsv);
    u(uid) = u(uid) - app.dt*dqdx;
    v(vid) = v(vid) - app.dt*dqdy;
    
    %correct pressure
    P(pid) = P(pid) + q(pid) - app.nu*divU;
    %remove mean pressure
    if imposemean
        P(pid)=P(pid)-mean(P(pid));
    else
        P([1:length(bcsDP),pid])=...
            P([1:length(bcsDP),pid])-mean(P(pid));
    end
    
    %Correct pressure for ui, vi's
    for i=1:app.S
        %Calculate RHS vector for P
        divU = Div2D(ui(:,i),vi(:,i),Nodeu,Nodev,dx,dy,idsP,Neu,Nev);
        RHSq = rhsPbc - 1/app.dt*divU;
        q(pid(Pperm))=UP\(LP\RHSq(Pperm));
        
        %Correct Velocity
        [dqdx dqdy]   = Grad2D(q,NodeP,dx,dy,idsu,idsv);
        ui(uid,i) = ui(uid,i) - app.dt*dqdx;
        vi(vid,i) = vi(vid,i) - app.dt*dqdy;
        
        %correct pressure
        Pi(pid,i) = Pi(pid,i) + q(pid) - app.nu*divU;
        %remove mean pressure
        if imposemean
            Pi(pid,i)=Pi(pid,i)-mean(Pi(pid,i));
        else
            Pi([1:length(bcsDP),pid],i)=...
                Pi([1:length(bcsDP),pid],i)-mean(Pi(pid,i));
        end
    end
    
    %%% CORRECT DO CONDITION FOR IMPLICIT TERMS
    %%%The number of iterations
    for iter = 1:app.DO_cor_iter
        %Calculate the correction for each ui, vi, and rhoi
        %use fui, fvi, and frhoi to save some memory
        for i = 1:app.S
            f2u = AU * ui00(uid, i);
            f2v = AV * vi00(vid, i);
            %Calculate inner product terms <ffiUj> -- there has to be app.S of them
            InProdi = zeros(app.S(1),1);
            for j = 1:app.S(1)
                InProdi(j) = (f2u' * ui00(uid, j) + ...
                    f2v' * vi00(vid, j)) * dx * dy;
            end
            fui(:, i) = ui00(uid, :) * InProdi;
            fvi(:, i) = vi00(vid, :) * InProdi;

            %Now repeat at the new time-level
            f2u = AU * ui(uid, i);
            f2v = AV * vi(vid, i);
            %Calculate inner product terms <ffiUj> -- there has to be app.S of them
            InProdi = zeros(app.S(1),1);
            for j = 1:app.S(1)
                InProdi(j) = (f2u' * ui(uid, j) + ...
                    f2v' * vi(vid, j)) * dx * dy;
            end
            fui(:, i) = fui(:, i) - ui(uid, :) * InProdi;
            fvi(:, i) = fvi(:, i) - vi(vid, :) * InProdi;
        end

        %Now do the corrections
        for i = 1:app.S
            ui(uid, i) = ui(uid, i) + app.dt * fui(:, i);
            vi(vid, i) = vi(vid, i) + app.dt * fvi(:, i);
        end
        %Record old basis for next iteration
        if iter < app.DO_cor_iter
            ui00 = ui;
            vi00 = vi;
        end
        %done, repeat if necessary
    end
    
    %% Ortho-Normalization of basis fields to avoid round-off errors
    [YYt ui vi Pi] = DOorthnorm(app, dx, dy, YYt, ui, uid, vi, vid, Pi, pid);
 
    %Advect Tracers (re-use RHSq to save a little memory)
    RHSq=1/app.dt*rho(pid) + app.kappa*rhsrho +...
        SCAadvect_TVD_mex(app.Nx, app.Ny, dx, dy, app.dt, rho, u, v, NodeP, Nodeu, Nodev, idsP,Nbcs,length(bcsDrho),NbcsDu,NbcsDv);
    %          SCAadvect_mex(app.Nx, app.Ny, dx, dy, rho, u, v, NodeP, Nodeu, Nodev, idsP,Nbcs,length(bcsDrho),NbcsDu,NbcsDv);
    rho(pid(rhoperm)) = Urho\(Lrho\RHSq(rhoperm));
    
    %Adapt stochastic dimension if required
    if ~mod(k,app.DOAdapt)
        DOAdapt_mex;
    end
    
    %Plot solution
    if ~mod(k,app.PlotIntrvl(1))
        eval(PlotScript);drawnow
    end
end
if mod(k,app.PlotIntrvl(1))
    eval(PlotScript);drawnow
end
%Clean up any temporary files
if exist('existbasis','var')
    delete([d0,sprintf('/Save/DOAdaptSave/Basis/orth_basis%0.3d.mat',existbasis)]);
end
