function [P u v rho NodeP Nodeu Nodev] =...
    FluidsSolver(app, SetupScript, PlotScript)
%function [P u v rho NodeP Nodeu Nodev] =...
%    FluidsSolver(app, SetupScript, PlotScript)
% Solves the Boussinesq(default)/Navier-Stokes/Stokes equations using a 
% simple Finite Volume discretization with an Incremental
% Pressure Correction in Rotational form projection method (by default).
% This function uses a second-order backward difference integration scheme,
% and uses LU factorization to solve the matrices. The advection uses a TVD
% advection scheme (default), and the advection is 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.PlotIntrvl: Number of time integration intervals between of 
%                          calling "PlotScript."
%          app.Stokes:     If true, this function solves the stokes
%                           equation
%          app.NonIncrement: If true, solves the equations using the 
%                           non-incremental Projection Method
%          app.NonRotational: if true, solves the equations using the
%                           Incremental Projection method in the 
%                           non-rotational form.
%          app.NonBoussinesq: If true, solves the Navier-Stokes instead of
%                           the Boussinesq equations
%          app.Advect: 'TVD' -- Use the Total Variation Diminishing 
%                               advection scheme (DEFAULT)
%                      'QUICK' -- Quadratic Upwind Interpolation scheme for
%                                   for Convective Kinetics
%                      'UW'  -- Use the UpWind advection scheme
%                      'CDS' -- Uses the Central Differencing Scheme
%           app.Fcor: Function that sets the coriolis forces. 
%               Called in the codes as: app.Fcor(time, [X, Y])
%       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
%           updateBcs: Function that updates the boundary conditions if
%                      these are time-dependent (optional). The Function
%                      call in the codes is as follows:
%            [bcsDP, bcsNP, bcsDrho, bcsNrho, bcsDu, bcsNu,bcsDv, bcsNv]...
%              = updateBcs(T, XP, YP, XU, YU, XV, YV, app, P, rho, u, v,...
%              NodeP, NodeRho, Nodeu, Nodev);
%           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

%Run setup script
eval(SetupScript)

%Options related to the Boussinesq equations
if ~exist('bcsDrho', 'var'), bcsDrho = zeros(size(bcsDP));end
if ~exist('bcsNrho', 'var'), bcsNrho = zeros(size(bcsNP));end
if ~exist('NodeRho', 'var'), NodeRho = NodeP;end
if ~exist('rho', 'var'), rho = zeros(size(P));end
if ~isfield(app, 'kappa'), app.kappa = 0;end
if ~exist('NodeRho', 'var'), NodeRho = NodeP;end
if ~isfield(app,'PrRa')
    app.PrRa=1;
end

%Options related to the Coriolis forces
if ~isfield(app,'Fcor'), app.Fcor = @(t, p) zeros(size(p, 1), 1);end

%Options related to boundary conditions
if ~exist('bcsOP','var'), bcsOP = [];end
if ~exist('bcsOu','var'), bcsOu = [];end
if ~exist('bcsOv','var'), bcsOv = [];end

%Options related to the solver
if ~isfield(app, 'Stokes'), app.Stokes = 0;end
if ~isfield(app, 'NonIncrement'), app.NonIncrement = 0;end
if ~isfield(app, 'NonRotational'), app.NonRotational = 0;end
if ~isfield(app, 'NonBoussinesq'), app.NonBoussinesq = 0;end
SlvOpts = [app.Stokes, app.NonIncrement,...
    app.NonRotational, app.NonBoussinesq];

%Options related to advection
if ~isfield(app, 'Advect'), app.Advect = 'TVD'; end

%Options related to boundary conditions
if exist('updateBcs','var') || exist('updateBcs','file')
    updatebcs=true;
else
    updatebcs=false;
end

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


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

[AP rhsPbc RHSP] = FixDiffPObcs(AP, bcsOP, Nbcs,...
    rhsPbc, NodeP, dx, dy, RHSP); bcsDP(bcsOP) = [];
[Au rhsubc RHSubc] = FixDiffPObcs(Au, bcsOu, Nbcs,...
    rhsubc, Nodeu, dx, dy, RHSubc); bcsDu(bcsOu) = [];
[Av rhsvbc RHSvbc] = FixDiffPObcs(Av, bcsOv, Nbcs,...
    rhsvbc, Nodev, dx, dy, RHSvbc); bcsDv(bcsOv) = [];

[Au Av rhsubc rhsvbc RHSubc RHSvbc] = FixDiffuvDbcs(Au, Av, bcsDu, bcsDv,...
    rhsubc, rhsvbc, Nodeu, Nodev, dx, dy, Nbcs, RHSubc, RHSvbc);

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

%And construc the Coriolis operators
if (sum(abs(app.Fcor(0, [XV(idsv) YV(idsv)])))>0)
    [Au2v Av2u] = CoriolisAVG(app.Nx, app.Ny, Nodeu, Nodev,...
        idsu, idsv, Nbcs);
    Au2v = spdiags(app.Fcor(0, [XV(idsv) YV(idsv)]), 0, ...
        length(idsv), length(idsv)) * Au2v;
    Av2u = spdiags(app.Fcor(0, [XU(idsu) YU(idsu)]), 0, ...
        length(idsu), length(idsu)) * Av2u;
else
    Au2v = spalloc(length(v) - Nbcs,length(u),0);
    Av2u = spalloc(length(u) - Nbcs,length(v),0);
end
% [Au2v Av2u] = CoriolisAVG(app.Nx, app.Ny, Nodeu, Nodev, idsu, idsv, Nbcs);
% 
% Au2v = spdiags(app.Fcor(0, [XV(idsv) YV(idsv)]), 0, length(idsv),...
%     length(idsv)) * Au2v;
% Av2u = spdiags(app.Fcor(0, [XU(idsu) YU(idsu)]), 0, length(idsu),...
%     length(idsu)) * Av2u;

%% 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);

% % LU factorization
[Lrho Urho] = lu(Arho(rhoperm,rhoperm));
[LP UP] = lu(AP(Pperm,Pperm));
[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];
NbcsDu = length(bcsDu); %Needed for advection
NbcsDv = length(bcsDv); %Needed for advection

%% Initialize Helper arrays
fu = zeros(length(uid), 1);
fv = zeros(length(vid), 1);

%% Solving system
q = zeros(size(P));
T = 0; k = 0; eval(PlotScript); drawnow
display('Paused at initial condition. Hit any key to continue');
drawnow; pause;

%Start with Implict Euler step
for k = 1:1
    T = app.dt * k;
    %Set boundary conditions
    if updatebcs
        [bcsDP, bcsNP, bcsDrho, bcsNrho, bcsDu, bcsNu,bcsDv, bcsNv]...
            = updateBcs(T, XP, YP, XU, YU, XV, YV, app, P, rho, u, v,...
            NodeP, NodeRho, Nodeu, Nodev);
        rho(1:Nbcs, :) = [bcsDrho'; bcsNrho'];
        rhsrho = -RHSrho * rho(1:Nbcs, :);
        P(1:Nbcs, :) = [bcsDP'; bcsNP'];
        rhsPbc = -RHSP * P(1:Nbcs, :);
        u(1:Nbcs, :) = [bcsDu'; bcsNu'];
        rhsubc = -RHSubc * u(1:Nbcs, :);
        v(1:Nbcs, :) = [bcsDv'; bcsNv'];
        rhsvbc = -RHSvbc * v(1:Nbcs, :);
    end    
    
    %If non-Incremental, do the next line
    if ~SlvOpts(2)
        [dpdx dpdy]   = Grad2D(P, NodeP, dx, dy, idsu, idsv);
    elseif k == 1
        dpdx = 0;
        dpdy = 0;
    else 
        dpdx = 0 * dpdx;
        dpdy = 0 * dpdy;
    end
    
    %If not the stokes solver, do advection
    if ~SlvOpts(1)
    switch lower(app.Advect)
        case 'tvd'
            [fu fv] = UVadvect_TVD_mex(app.Nx, app.Ny, dx, dy, app.dt,...
                u, v, u, v, Nodeu, Nodev, idsu, idsv, Nbcs, NbcsDu, NbcsDv);
        case 'quick'
            [fu fv] = UVadvect_DO_mex(app.Nx, app.Ny, dx, dy, ...
                u, v, u, v, Nodeu, Nodev, idsu, idsv, Nbcs, NbcsDu, NbcsDv);
        case 'uw'
            [fu fv] = UVadvect_UW(app.Nx, app.Ny, dx, dy, ...
                u, v, Nodeu, Nodev, idsu, idsv);
        case 'cds'
            [fu fv] = UVadvect_DO_CDS_mex(app.Nx, app.Ny, dx, dy,...
            u, v, u, v, Nodeu, Nodev, idsu, idsv, Nbcs, NbcsDu, NbcsDv);
    end
    end
    
    %If Boussinesq add the density to the source term
    if ~SlvOpts(4)
        src = app.PrRa * 0.5 * ( rho(NodeRho(srcTopid, srcid))...
            + rho(NodeRho(srcBotid, srcid)) );
        src = src(idsv);
    else 
        src = 0;
    end
    
    RHSu = 1 / app.dt * u(uid) + fu + rhsubc * app.nu - dpdx + ...
        Fu(T, [XU(idsu) YU(idsu)]) + Av2u*v;
    RHSv = 1 / app.dt * v(vid) + fv + rhsvbc * app.nu - dpdy - src + ...
        Fv(T, [XV(idsv) YV(idsv)]) - Au2v*u;

    %Save old values
    uold = u;
    vold = v;

    %For this timestep, we need to remove 1/2 app.dt from the diagonal of
    % Au, Av for the first order Euler timestepping
    u(uid) = (Au - 0.5 / app.dt * speye(size(Au))) \ RHSu;
    v(vid) = (Av - 0.5 / app.dt * speye(size(Av))) \ RHSv;

    %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
    [dpdx dpdy]   = Grad2D(q, NodeP, dx, dy, idsu, idsv);
    u(uid) = u(uid) - app.dt * dpdx;
    v(vid) = v(vid) - app.dt * dpdy;

    %correct pressure
    if ~SlvOpts(2)
        P(pid) = P(pid) + q(pid) - app.nu * divU * (1  -SlvOpts(3));
    end
    
    %remove mean pressure
    if imposemean
        P = P - mean(P(:));
    end
    
    %Advect Tracers
    switch lower(app.Advect)
        case 'tvd'
            RHSq = SCAadvect_TVD_mex(app.Nx, app.Ny, dx, dy, app.dt,...
                rho, u, v, NodeRho, Nodeu, Nodev, idsP, Nbcs,...
                length(bcsDrho), NbcsDu, NbcsDv);
        case 'quick'
            RHSq = SCAadvect_mex(app.Nx, app.Ny, dx, dy,...
                rho, u, v, NodeRho, Nodeu, Nodev, idsP, Nbcs, ...
                length(bcsDrho), NbcsDu, NbcsDv);
        case 'uw'
            RHSq = SCAadvect_UW(app.Nx, app.Ny, dx, dy, rho, u, v,...
                NodeRho, Nodeu, Nodev, idsP);
        case 'cds'
            RHSq = SCAadvect_CDS_mex(app.Nx, app.Ny, dx, dy, ...
                rho, u, v, NodeP, Nodeu, Nodev, idsP, Nbcs, ...
                length(bcsDrho), NbcsDu, NbcsDv);
    end
    RHSq = RHSq + 1 / app.dt * rho(pid) + app.kappa * rhsrho;
    %Save old values
    rhoold = rho;
    rho(pid) = (Arho - 0.5 / app.dt * speye(size(Arho))) \ RHSq;
    
    %Plot solution
    if ~mod(k,app.PlotIntrvl)
        eval(PlotScript);drawnow
    end
end

%Continue with second order BDF scheme
for k = 2:Nt
    T = app.dt * k;
    
    %Set boundary conditions
    if updatebcs
        [bcsDP, bcsNP, bcsDrho, bcsNrho, bcsDu, bcsNu,bcsDv, bcsNv]...
            = updateBcs(T, XP, YP, XU, YU, XV, YV, app, P, rho, u, v,...
            NodeP, NodeRho, Nodeu, Nodev);
        rho(1:Nbcs, :) = [bcsDrho'; bcsNrho'];
        rhsrho = -RHSrho * rho(1:Nbcs, :);
        P(1:Nbcs, :) = [bcsDP'; bcsNP'];
        rhsPbc = -RHSP * P(1:Nbcs, :);
        u(1:Nbcs, :) = [bcsDu'; bcsNu'];
        rhsubc = -RHSubc * u(1:Nbcs, :);
        v(1:Nbcs, :) = [bcsDv'; bcsNv'];
        rhsvbc = -RHSvbc * v(1:Nbcs, :);
    end
            
    %If non-Incremental, do the next line
    if ~SlvOpts(2)
        [dpdx dpdy]   = Grad2D(P, NodeP, dx, dy, idsu, idsv);
    elseif k == 1
        dpdx = 0;
        dpdy = 0;
    else 
        dpdx = 0 * dpdx;
        dpdy = 0 * dpdy;
    end
    
    %If not the stokes solver, do advection
    if ~SlvOpts(1)
    switch lower(app.Advect)
        case 'tvd'
            [fu fv] = UVadvect_TVD_mex(app.Nx, app.Ny, dx, dy, app.dt,...
                u, v, u, v, Nodeu, Nodev, idsu, idsv, Nbcs, NbcsDu, NbcsDv);
        case 'quick'
            [fu fv] = UVadvect_DO_mex(app.Nx, app.Ny, dx, dy, ...
                u, v, u, v, Nodeu, Nodev, idsu, idsv, Nbcs, NbcsDu, NbcsDv);
        case 'uw'
            [fu fv] = UVadvect_UW(app.Nx, app.Ny, dx, dy, ...
                u, v, Nodeu, Nodev, idsu, idsv);
        case 'cds'
            [fu fv] = UVadvect_DO_CDS_mex(app.Nx, app.Ny, dx, dy,...
                u, v, u, v, Nodeu, Nodev, idsu, idsv, Nbcs, NbcsDu, NbcsDv);
    end
    end

    %If Boussinesq add the density to the source term
    if ~SlvOpts(4)
        src = app.PrRa * 0.5 * ( rho(NodeRho(srcTopid, srcid))...
            + rho(NodeRho(srcBotid, srcid)) );
        src = src(idsv);
    else 
        src = 0;
    end    
    
    RHSu = 2 / app.dt * u(uid) - 1 / (2 * app.dt) * uold(uid) + fu +...
        rhsubc * app.nu - dpdx +...
        Fu(T, [XU(idsu) YU(idsu)]) + Av2u*v;
    RHSv = 2 / app.dt * v(vid) - 1 / (2 * app.dt) * vold(vid) + fv + ...
        rhsvbc * app.nu - dpdy - src + ...
        Fv(T, [XV(idsv) YV(idsv)]) - Au2v*u;
    
    %Save old values
    uold = u;
    vold = v;
    
    u(uid(uperm)) = Uu \ (Lu \ RHSu(uperm));
    v(vid(vperm)) = Uv \ (Lv \ RHSv(vperm));

    %Calculate RHS vector for P
    divU = Div2D(u, v, Nodeu, Nodev,dx, dy, idsP, Neu, Nev);
    RHSq = rhsPbc - 3 / 2 / app.dt * divU;
    q(pid(Pperm)) = UP \ (LP \ RHSq(Pperm));

    %Correct Velocity
    [dpdx dpdy]   = Grad2D(q, NodeP, dx, dy, idsu, idsv);
    u(uid) = u(uid) - 2 / 3 * app.dt * dpdx;
    v(vid) = v(vid) - 2 / 3 * app.dt * dpdy;

    %correct pressure
    if ~SlvOpts(2)
        P(pid) = P(pid) + q(pid) - app.nu * divU * (1 - SlvOpts(3));
    end
    
    %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
    
    %Advect Tracers
    switch lower(app.Advect)
        case 'tvd'
            RHSq = SCAadvect_TVD_mex(app.Nx, app.Ny, dx, dy, app.dt,...
                rho, u, v, NodeRho, Nodeu, Nodev, idsP, Nbcs,...
                length(bcsDrho), NbcsDu, NbcsDv);
        case 'quick'
            RHSq = SCAadvect_mex(app.Nx, app.Ny, dx, dy,...
                rho, u, v, NodeRho, Nodeu, Nodev, idsP, Nbcs, ...
                length(bcsDrho), NbcsDu, NbcsDv);
        case 'uw'
            RHSq = SCAadvect_UW(app.Nx, app.Ny, dx, dy, rho, u, v,...
                NodeRho, Nodeu, Nodev, idsP);
        case 'cds'
            RHSq = SCAadvect_CDS_mex(app.Nx, app.Ny, dx, dy, ...
                rho, u, v, NodeP, Nodeu, Nodev, idsP, Nbcs, ...
                length(bcsDrho), NbcsDu, NbcsDv);
    end
    RHSq = RHSq + 2 / app.dt * rho(pid) - 1 / (2 * app.dt) * rhoold(pid)...
        + app.kappa * rhsrho;
    
    %Save old values
    rhoold = rho;
    rho(pid(rhoperm)) = Urho \ (Lrho \ RHSq(rhoperm));

    %Plot solution
    if ~mod(k,app.PlotIntrvl)
        eval(PlotScript);drawnow
    end
end
if mod(k, app.PlotIntrvl)
    eval(PlotScript); drawnow
end
