function [P, u, v, rho, NodeP, Nodeu, Nodev] = FluidsSolver2(app, SetupScript, PlotScript)
%function [P u v rho NodeP Nodeu Nodev] = FluidsSolver2(app, SetupScript, PlotScript)
%
% Solves the Navier-Stokes/Stokes equations with Boussinesq buoyancy approximation by
% a simple Finite Volume discretization on a uniform Cartesian mesh using
% an Incremental Pressure Correction in Rotational Form projection method (by default).
% This function uses a second-order backward difference time marching scheme,
% and uses LU factorization to solve the linear systems. The advection uses a TVD
% advection scheme (default), and the advection is mexed.
%
%
% INPUTS:
%
% app:	          The application data structure of which the fields can also be
%                 added through the SetupScript.
% SetupScript:    String input of filename used to setup the problem.
% 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.
%
%
% At a minimum, the structure app should contain the following fields:
%
% app.Nx:         Number of interior P cells in x-direction (integer number)
% app.Ny:         Number of interior P cells in y-direction (integer number)
% app.T:          Final time
% app.dt:         Timestep size
% app.PlotIntrvl: Number of time steps between two "PlotScript" calls
%
% Other optional fields of app:
%
% app.nu:         Effective viscosity in x direction
% app.nu2:        Effective viscosity in y direction
% app.kappa:      Effective diffusivity of density in x direction
% app.kappa2:     Effective diffusivity of density in y direction
% app.g:          Effective gravity acceleration along negative y direction
% app.Stokes:     If true, this function solves the Stokes equation
% app.NonInc:     If true, solves the equations using the
%                 non-incremental Projection Method
% app.NonRot:     If true, solves the equations using the
%                 Incremental Projection method in the non-rotational form
% 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'   -- Use the Central Differencing Scheme
% app.Fcor:       Function that sets the coriolis forces.
%                 Called in the codes as: app.Fcor(time, [X, Y])
% app.updateBcs:  (Optional) Function that updates the time dependent boundary conditions.
%
%
% This solver needs the following variables correctly created in the SetupScript.
%
% 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
% bcsDq:     q (for Projection Method) Dirichlet boundary conditions
% bcsDu:     U-velocity Dirichlet boundary conditions
% bcsDv:     V-velocity Dirichlet boundary conditions
% bcsNrho:   Density Neumann boundary conditions
% bcsNq:     q (for Projection Method) 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 cell id's (see NodePad.m)
% idsu:      Interior U-velocity cell id's (see NodePad.m)
% idsv:      Interior V-velocity cell id's (see NodePad.m)
% Fu(t, XY): Forcing function for U-velocity, XY =[XU(idsu), YU(idsu)]
% Fv(t, XY): Forcing function for V-velocity, XY =[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 for q.
% 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 for q.
% 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: Diffusion_Operator, Grad2D, Div2D
%
% Date: Jan. 2014
% Author: Matt Ueckermann, Jing Lin

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Preliminary setup

% Run setup script
eval(SetupScript);

% Default options related to the solver
if ~isfield(app, 'Stokes'), app.Stokes = 0; end;
if ~isfield(app, 'NonInc'), app.NonInc = 0; end;
if ~isfield(app, 'NonRot'), app.NonRot = 0; end;
SlvOpts = [app.Stokes, app.NonInc, app.NonRot];

% Default options related to the tracer equations
if ~exist('bcsDrho', 'var'), bcsDrho = zeros(size([bcsDq; bcsNq])); end;
if ~exist('bcsOrho', 'var'), bcsOrho = []; end;
if ~exist('bcsNrho', 'var'), bcsNrho = []; 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 ~isfield(app, 'kappa2'), app.kappa2 = app.kappa; end;
if ~isfield(app,'g'), app.g = 1; end;
if ~isfield(app,'Frho'), app.Frho = @(t, XY, rho) zeros(size(XY, 1), 1); end;

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

% Default options related to diffusion
if ~isfield(app, 'nu'), app.nu = 1; end;
if ~isfield(app, 'nu2'), app.nu2 = app.nu; end;

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

% Default options related to open boundary conditions
if ~exist('bcsOq','var'), bcsOq = []; end;
if ~exist('bcsOu','var'), bcsOu = []; end;
if ~exist('bcsOv','var'), bcsOv = []; end;
if ~exist('bcsOrho','var'), bcsOrho = []; end;
if (~isempty(bcsOu))||(~isempty(bcsOv))||(~isempty(bcsOrho))
    error('Open boundary conditions can only be set for pressure!');
end;

% Find number of boundary conditions and number of interior cells
Nbcs = length(bcsDq) + length(bcsNq);
Num_IntCell_P = max(NodeP(:)) - Nbcs;
Num_IntCell_u = max(Nodeu(:)) - Nbcs;
Num_IntCell_v = max(Nodev(:)) - Nbcs;


%% Construct diffusion operator matrices and corresponding RHS vectors

[Dfs_rho, Dfs_bc_rho, Dfs_bc_rhs_rho] = Diffusion_Operator...
    ('P', NodeRho, bcsDrho, bcsOrho, bcsNrho, dx, dy, app.kappa, app.kappa2);
[Dfs_q, Dfs_bc_q, Dfs_bc_rhs_q] = Diffusion_Operator...
    ('P', NodeP, bcsDq, bcsOq, bcsNq, dx, dy, 1, 1, true);
[Dfs_u, Dfs_bc_u, Dfs_bc_rhs_u] = Diffusion_Operator...
    ('u', Nodeu, bcsDu, bcsOu, bcsNu, dx, dy, app.nu, app.nu2);
[Dfs_v, Dfs_bc_v, Dfs_bc_rhs_v] = Diffusion_Operator...
    ('v', Nodev, bcsDv, bcsOv, bcsNv, dx, dy, app.nu, app.nu2);

Dfs_rho = Dfs_rho + 3/2/app.dt * speye(size(Dfs_rho));
Dfs_u = Dfs_u + 3/2/app.dt * speye(size(Dfs_u));
Dfs_v = Dfs_v + 3/2/app.dt * speye(size(Dfs_v));

% Construct 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


%% LU factorization of matrices

% First order matrix to give best spartsity pattern
rhoperm = symamd(Dfs_rho);
qperm = symamd(Dfs_q);
uperm = symamd(Dfs_u);
vperm = symamd(Dfs_v);

% LU factorization
[Lrho, Urho] = lu(Dfs_rho(rhoperm,rhoperm));
[Lq, Uq] = lu(Dfs_q(qperm,qperm));
[Lu, Uu] = lu(Dfs_u(uperm,uperm));
[Lv, Uv] = lu(Dfs_v(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;
bcN_ID_u = [length(bcsDu)+1, Nbcs];
bcN_ID_v = [length(bcsDv)+1, Nbcs];
NbcsDu = length(bcsDu); % Needed for advection
NbcsDv = length(bcsDv); % Needed for advection


%% Time marching to solve the time-evolving system

tic;

% Helper arrays which store advection terms
fu = zeros(length(uid), 1);
fv = zeros(length(vid), 1);

q = zeros(size(P));
T = 0; k = 0;

% Remove the mean of the pressure (determined up to a constant)
P(pid) = P(pid) - mean(P(pid));

eval(PlotScript);
fprintf('Paused at initial condition. Hit any key to continue.\n\n');
pause;

% Use backward Euler (BD1) scheme for the first time step

for k = 1 : 1
    T = app.dt * k;
    
    % Set time-dependent boundary conditions
    if isfield(app, 'updateBcs')
        [bcsDrho, bcsNrho, bcsDu, bcsNu, bcsDv, bcsNv]...
            = app.updateBcs(T, XP, YP, XU, YU, XV, YV, app, NodeRho, Nodeu, Nodev,...
            bcsDrho, bcsNrho, bcsDu, bcsNu, bcsDv, bcsNv);
        rho(1:Nbcs) = [bcsDrho; bcsNrho];
        Dfs_bc_rhs_rho = -Dfs_bc_rho * rho(1:Nbcs);
        u(1:Nbcs) = [bcsDu; bcsNu];
        Dfs_bc_rhs_u = -Dfs_bc_u * u(1:Nbcs);
        v(1:Nbcs) = [bcsDv; bcsNv];
        Dfs_bc_rhs_v = -Dfs_bc_v * v(1:Nbcs);
    end;
    
    % If the incremental projection method is used,
    % compute the gradient of the pressure in the last time step.
    if ~app.NonInc
        [dpdx, dpdy] = Grad2D_P(P, NodeP, dx, dy, idsu, idsv);
    else
        dpdx = 0; dpdy = 0;
    end;
    
    % If not the stokes solver, compute advection terms of momentum equations.
    if ~app.Stokes
        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;
    
    % Add the gravity source term
    src = app.g*0.5*(rho(NodeRho(srcTopid, srcid)) + rho(NodeRho(srcBotid, srcid)));
    src = src(idsv);
    
    RHSu = 1/app.dt * u(uid) + fu + Dfs_bc_rhs_u - dpdx...
        + Fu(T, [XU(idsu) YU(idsu)]) + Av2u*v;
    RHSv = 1/app.dt * v(vid) + fv + Dfs_bc_rhs_v - dpdy...
        + Fv(T, [XV(idsv) YV(idsv)]) - Au2v*u - src;
    
    % Save old values for BD2 time marching
    u_old = u; v_old = v;
    
    % Compute the intermediate velocity
    % For this timestep, we need to remove 1/(2*dt) from the diagonal entries of
    % Dfs_u, Dfs_v for the first time step using backward Euler(BD1).
    u(uid) = (Dfs_u - 0.5/app.dt*speye(size(Dfs_u))) \ RHSu;
    v(vid) = (Dfs_v - 0.5/app.dt*speye(size(Dfs_v))) \ RHSv;
    
    % Projection step (Helmholtz decomposition).
    % Compute the potential q of the curl-free part of the intermediate velocity
    [dudx, dvdy] = Grad2D_uv(u, v, Nodeu, Nodev, idsP, dx, dy, bcN_ID_u, bcN_ID_v);
    RHSq = Dfs_bc_rhs_q - 1/app.dt*(dudx + dvdy);
    q(pid(qperm)) = Uq \ (Lq \ RHSq(qperm));
    
    % Correct Velocity
    [dpdx, dpdy] = Grad2D_P(q, NodeP, dx, dy, idsu, idsv);
    u(uid) = u(uid) - app.dt * dpdx;
    v(vid) = v(vid) - app.dt * dpdy;
    
    % Correct pressure
    if ~app.NonInc
        P(pid) = P(pid) + q(pid);
        if ~app.NonRot
            P(pid) = P(pid) - (app.nu*dudx + app.nu2*dvdy);
            [dudx, dvdy] = Grad2D_uv(u, v, Nodeu, Nodev, idsP, dx, dy, bcN_ID_u, bcN_ID_v);
            P(pid) = P(pid) + (app.nu*dudx + app.nu2*dvdy);
        end;
    end;
    
    % Remove the mean of the pressure (determined up to a constant)
    P(pid) = P(pid) - mean(P(pid));
    
    % Compute tracer advection terms
    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, NodeRho, Nodeu, Nodev, idsP, Nbcs, ...
                length(bcsDrho), NbcsDu, NbcsDv);
    end;
    
    RHSq = RHSq + 1/app.dt * rho(pid) + Dfs_bc_rhs_rho...
        + app.Frho(T, [XP(idsP) YP(idsP)], rho(pid));
    
    % Save old values for BD2 time marching
    rho_old = rho;
    rho(pid) = (Dfs_rho - 0.5/app.dt * speye(size(Dfs_rho))) \ RHSq;
    
    % Plot and save data to file
    if ~mod(k,app.PlotIntrvl), eval(PlotScript); end;
end;

% Use second order backward differencing (BD2) scheme from the second time step

for k = 2 : Nt
    T = app.dt * k;
    
    % Set time-dependent boundary conditions
    if isfield(app, 'updateBcs')
        [bcsDrho, bcsNrho, bcsDu, bcsNu, bcsDv, bcsNv]...
            = app.updateBcs(T, XP, YP, XU, YU, XV, YV, app, NodeRho, Nodeu, Nodev,...
            bcsDrho, bcsNrho, bcsDu, bcsNu, bcsDv, bcsNv);
        rho(1:Nbcs) = [bcsDrho; bcsNrho];
        Dfs_bc_rhs_rho = -Dfs_bc_rho * rho(1:Nbcs);
        u(1:Nbcs) = [bcsDu; bcsNu];
        Dfs_bc_rhs_u = -Dfs_bc_u * u(1:Nbcs);
        v(1:Nbcs) = [bcsDv; bcsNv];
        Dfs_bc_rhs_v = -Dfs_bc_v * v(1:Nbcs);
    end;
    
    % If the incremental projection method is used,
    % compute the gradient of the pressure in the last time step.
    if ~app.NonInc
        [dpdx, dpdy] = Grad2D_P(P, NodeP, dx, dy, idsu, idsv);
    else
        dpdx = 0; dpdy = 0;
    end;
    
    % If not the stokes solver, compute advection terms of momentum equations.
    if ~app.Stokes
        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;
    
    % Add the gravity source term
    src = app.g*0.5*(rho(NodeRho(srcTopid, srcid)) + rho(NodeRho(srcBotid, srcid)));
    src = src(idsv);
    
    RHSu = 2/app.dt * u(uid) - 1/(2*app.dt) * u_old(uid) + fu + Dfs_bc_rhs_u - dpdx...
        + Fu(T, [XU(idsu) YU(idsu)]) + Av2u*v;
    RHSv = 2/app.dt * v(vid) - 1/(2*app.dt) * v_old(vid) + fv + Dfs_bc_rhs_v - dpdy...
        + Fv(T, [XV(idsv) YV(idsv)]) - Au2v*u - src;
    
    % Save old values for BD2 time marching
    u_old = u; v_old = v;
    
    % Compute the intermediate velocity
    u(uid(uperm)) = Uu \ (Lu \ RHSu(uperm));
    v(vid(vperm)) = Uv \ (Lv \ RHSv(vperm));
    
    % Projection step (Helmholtz decomposition).
    % Compute the potential q of the curl-free part of the intermediate velocity
    [dudx, dvdy] = Grad2D_uv(u, v, Nodeu, Nodev, idsP, dx, dy, bcN_ID_u, bcN_ID_v);
    RHSq = Dfs_bc_rhs_q - 3/2/app.dt*(dudx + dvdy);
    q(pid(qperm)) = Uq \ (Lq \ RHSq(qperm));
    
    % Correct velocity
    % Subtract the curl-free velocity from the intermedate one to obtain the desired
    % divergence-free part.
    [dpdx, dpdy]   = Grad2D_P(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 ~app.NonInc
        P(pid) = P(pid) + q(pid);
        if ~app.NonRot
            P(pid) = P(pid) - (app.nu*dudx + app.nu2*dvdy);
            [dudx, dvdy] = Grad2D_uv(u, v, Nodeu, Nodev, idsP, dx, dy, bcN_ID_u, bcN_ID_v);
            P(pid) = P(pid) + (app.nu*dudx + app.nu2*dvdy);
        end;
    end;
    
    % Remove the mean of the pressure (determined up to a constant)
    P(pid) = P(pid) - mean(P(pid));
    
    % Compute tracer advection terms
    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, NodeRho, Nodeu, Nodev, idsP, Nbcs, ...
                length(bcsDrho), NbcsDu, NbcsDv);
    end;
    RHSq = RHSq + 2/app.dt*rho(pid) - 1/(2*app.dt)*rho_old(pid) + Dfs_bc_rhs_rho...
        + app.Frho(T, [XP(idsP) YP(idsP)], rho(pid));
    
    % Save old values for BD2 time marching
    rho_old = rho;
    rho(pid(rhoperm)) = Urho \ (Lrho \ RHSq(rhoperm));
    
    % Plot and save data to file
    if ~mod(k,app.PlotIntrvl), eval(PlotScript); end;
end;

