%% Scalar advection demo (1st order Upwind)

% Last modified by Jing Lin in Apr. 2014


%% Set up data and images saving directory


%% Set up problem-specific parameterr.


%% Set up mesh and solution domain

% Define size of solution domain
a = 1;
b = 1;

if ~isfield(app, 'Nx'), error('app.Nx not specified'); end;
if ~isfield(app, 'Ny'), error('app.Ny not specified'); end;

% Set number of timesteps.
if ~isfield(app, 'T'), error('app.T not specified'); end;
if ~isfield(app, 'dt'), error('app.dt not specified'); end;
Nt = ceil(app.T/app.dt);
if ~isfield(app, 'PlotIntrvl'), error('app.PlotIntrvl not specified'); end;

% Mesh the domain
% This is a standard call and should not require much modification.
% [0, a] gives the min/max x coordinates
% and [0, b] gives the min/max y coordinates.
[XP, YP, XU, YU, XV, YV, dx, dy] = meshdomain(app, [0, a], [0, b]);


%% Set up the geometry mask

% This mask is user-created. Zeros are unmasked, and non-zero integers
% indicate a masked region. The domain sides should not be considered here.
% Each interior masked region could have a different number to allow
% different boundary conditions.
% Create mask of correct size. This is a standard call and should not
% require modification.
Mask = zeros(app.Ny, app.Nx);

% Logical array add_bnd signifies whether to add (add if =true) boundaries 
% for the [west, east, south, north] sides of the solution domain.
% This is almost a standard procedure except when periodic boundaries are used.
add_bnd = [false, false, false, false];

% Set periodic BC flag
periodic = 1;

% Create preliminary node numbering matrices using function "NodePad".
% Note: This is a standard procedure and should not be modified.
[NodeP, Nodeu, Nodev, idsP, idsu, idsv] = NodePad(Mask, add_bnd, periodic);

% Plot the current boundary IDs to facilitate setting up boundary conditions.
if ~app.NoPlot
    plot_bc_ids(NodeP, idsP);
    fprintf('Paused at boundary id plot. Hit any key to continue.\n\n');
    pause;
end;

%% Set up boundary conditions (BCs)

% User assigns BC types to each boundary ID for u, v, rho in three cell arrays.

% This requires user input. Use the above plot as reference to identify BC IDs.
% e.g. Element 1 in cell array bcid2type_u represents the BC type of
% the boundary with ID 1 for velocity u.

% Note: Only the following four types of BCs can be set for velocities
% and density:
% Single-value Dirichlet BC   (='D')
% Multiple-value Dirichlet BC (='Dmv')
% Single-value Neumann BC     (='N')
% Multiple-value Neumann BC   (='Nmv')

% Warning: Open boundary conditions (='O') are not allowed for velocities
% and density!

bcid2type_u = {};      
bcid2type_v = {};
bcid2type_q = {};

% Set up the final node numbering matrix using function "set_bcs".
% For each single-value BC ('D', 'N', 'O'), all its boundary cells share one ID,
% like in the preliminary node numbering matrices.
% For each multiple-value BC ('Dmv', 'Nmv'), now each boundary cell will
% be assigned a different ID.
% Besides, the IDs are arranged to satisfy: 
%           Dirichlet_BC_IDs < Open_BC_IDs < Neumann_BC_IDs

% Note: This is a standard procedure and should not be modified.

% Without density equation:
[NodeP, nbcq, Nodeu, nbcu, Nodev, nbcv] = set_bcs(...
    bcid2type_q, NodeP, bcid2type_u, Nodeu, bcid2type_v, Nodev);

% With density equation:
% [NodeP, nbcq, Nodeu, nbcu, Nodev, nbcv, NodeRho, nbcrho] = set_bcs(... 
%     bcid2type_q, NodeP, bcid2type_u, Nodeu, bcid2type_v, Nodev, ...
%     bcid2type_rho, NodeP);

% Set preliminary BC values to zero since most BCs are homogeneous.
% Some of them can be modified to be non-homogeneous BCs later.

% Note: This is a standard procedure and should not be modified.

[bcsDq, bcsOq, bcsNq] = init_bcs(nbcq);
[bcsDu, bcsOu, bcsNu] = init_bcs(nbcu);
[bcsDv, bcsOv, bcsNv] = init_bcs(nbcv);

% User sets non-zero boundary conditions.


%% Set up initial conditions

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

P_init   = @(X, Y) sin(2*pi.*X/a).*sin(2*pi.*Y/b);
u_init   = @(X, Y) a.*ones(size(X));
v_init   = @(X, Y) b.*ones(size(X));

% Initialize the state variables with values from both boundary and internal
% cells.

% Note: This is a standard procedure and should not be modified.

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


%% Preliminary setup

% Find number of boundary conditions and number of interior cells
Nbcs = length(bcsDq) + length(bcsNq);

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


%% Plot the initial state

figure;
surf(XP, YP, P(NodeP(2:app.Ny+1,2:app.Nx+1)), 'LineStyle', 'none');
colorbar; zlim([-1 1]);
TT(1) = xlabel('x(m)'); TT(2) = ylabel('y(m)'); TT(3) = zlabel('T(^oC)');
TT(4) = title(sprintf('Numerical Solution Nx=%g, Ny=%g', app.Nx, app.Ny));
formatplot([],gca,TT); caxis([-1,1]);
drawnow;

display('Paused at initial condition. Hit any key to continue');
pause;


%% Time Marching

for k = 1 : Nt
    P(pid,1) = P(pid,1) + app.dt* ...
        SCAadvect_UW(app.Nx, app.Ny, dx, dy, P, u, v, NodeP, Nodeu, Nodev, idsP);
    if ~mod(k, app.PlotIntrvl)
        surf(XP, YP, P(NodeP(2:app.Ny+1,2:app.Nx+1)), 'LineStyle', 'none');
        colorbar; zlim([-1 1]);
        TT(1)=xlabel('x(m)'); TT(2)=ylabel('y(m)'); TT(3)=zlabel('T(^oC)');
        TT(4)=title(sprintf('Numerical Solution Nx=%g, Ny=%g',app.Nx,app.Ny));
        formatplot([],gca,TT); caxis([-1,1]);
        drawnow;
    end;
end;

